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Abstract 

DEVELOPMENT OF 3-D EM SIMULATOR FOR DESIGNING 
MICROWAVE CIRCUITS USING MATERIAL INDEPENDENT 
PERFECTLY MATCHED LAYERS AS ABC 

Thesis Supervisor Dr. Animesh Biswas 

Co-guide Dr. Joseph John 
Akhilesh Jain 


An improved version of unsplit perfectly matched layers (PML). which is material independent, has 
been formulated for the simulation of microwave circuits embedded in lossy, isotropic or anisotropic, 
dielectric and magnetic 3-D media. Advantage of using PML is that it provides virtually reflection 
free absorption, which gives highly accurate solutions for circuit parameters. It has been demonstrated 
that because of accuracy of material independent perfectly matched layers (MIPML), the usual 
oscillatory behavior observed in dispersion relations of characteristic impedance and relative effective 
dielectric constant in FDTD analysis is drastically reduced. The other advantage of MIPML 
formulation is that for various t>pes of geometries and material, there is no modification in the 
algorithm for the calculation of flux components in perfectly matched layers and in working domain. 
By adding a little structural details in the subroutines; responsible for setting up PEC in w'orking 
domain, for calculating E and H fields from flux components and for specifying source stimulus, the 
program can calculate circuit parameters. This greatly enhances potential for calculating responses of 
complex circuits and radiating structures without any major changes in the program. Though the 
scheme, which has been adapted here requires more memory but it lends it self for uniformity of the 
algorithm throughout the computing domain, reduces number of subroutines and makes it easier for ’ 
developers to debug the program. A collection of sample problems have been included to demonstrate 
the w'orking potential of the program and to validate calculated results with published data. Samples 
include patch antenna on an anisotropic substrate, propagation of a Gaussian pulse in a microstrip, 
wave propagating in a uniform anisotropic 3-D medium. It has been found that calculated and 
published results are in excellent agreement. Dispersion relationships of characteristic impedance and 
relative effective dielectric constants of microstrip on an arbitrary' anisotropic substrate, which are 
very useful for CAD have been worked out. Here these responses have been calculated by both, 
conventional method of multiple cells and a new single cell method based upon transmission line 
analogy. It has been successfully demonstrated that new method developed in this thesis further 
reduces oscillations observed in dispersion relations. Improvement is quite distinctive at lower 
frequency side. 



Chapter 1 
Introduction 

Finite difference time domain method was introduced by Yee K, S. [1] back in 1966. 
Method did not receive much attention due to several reasons. Prominent of them was 
high cost of computation. Although with time, cost of computation came down, there 
was a lack of proper absorbing boundary conditions, which were inherent in the original 
scheme for open boundary problems. After researchers proposed several ABCs, there has 
been a tremendous growth of activity in FDTD making it the most popular method in 
electromagnetics problems. There are numerous advantages of the FDTD method. In 
present scenario of computation, it enjoys popularity because it lends itself to be 
programmed on parallel processors and thus analysis of very complex geometries is 
possible. Further being a time domain technique, in a single run it gives response of the 
geometry over a wide frequency range. Another reason for its popularity is that very little 
analytical preparation is required and standard algorithms can be developed which cover 
a wide range of the problems. The scope of FDTD method covers almost entire field of 
electromagnetics, so it is difficult to cover all of them. The FDTD method has been 
extensively used for Microwave CAD, mostly for isotropic medium. In this thesis FDTD 
method has been chosen over other methods (TLM or Bergeron’s) because it is extremely 
efficient, its implementation is quite straight forward, and it may be derived directly from 
Maxwell’s equation. X.Zhang et al. [2,3] had used FDTD with Super Absorbing 
Boundary conditions to analyze frequency dependent characteristics of microstrip 
transmission lines and discontinuities. They had used electric wall source, which resulted 
in a sharp magnetic field tangential to the source wall. This results in some distortion of 
the launched pulse. Sheen et al [11], for characterizing microstrip discontinuities, also 
used FDTD method with Mur’s first order ABC and used magnetic wall in source plane, 
which was imposed using image theory which avoided induced magnetic field. Here 
electric field components on source plane were obtained using finite difference equations. 
However, method of excitation was complicated but it was better than Xhang et al. In 
addition, ABC used by them was simpler. In FDTD analysis, frequency dependence of 
microstrip parameters shows oscillatory behavior [12]. One of the contributors to 
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oscillatory behavior is mixed modes of propagation and evanescent waves contained in 
excitation source. Advanced excitation source scheme proposed by Zhao in [26,27] is 
quite simple and much better. Major contribution to oscillatory behavior comes from 
imperfect ABCs and the Mur’s ABC, used by above researchers is not perfect. Also 
computation domain required by these pioneering works was large as ABCs used, 
required scatterer to be placed quite far, so impinging waves on interface become plane 
waves and near normal incident. In October 1 994, Berenger created waves by introducing 
Perfectly Matched Layer [13], which was far superior, compared to existing ABCs. The 
innovation of Berenger’ s PML is that plane waves of arbitrary incidence, polarization, 
and frequency are matched at the boundary. This helped in reducing size of computing 
domain enormously. Though Berenger’ s original paper was for two-dimensional isotropic 
space but Katz [25] et al. validated it for 3-D space. Its wide band performance was 
validated by Reuter et al. [26]. Among other researchers who put forward rigorous 
analysis of absorption characteristics of PML, were Prescott et al. [29,30] who compared 
PML with other ABCs in a two-part article. Extensive research by Fang et al. [15,16] 
overcame certain disadvantages of PML formulation, mainly absorption of evanescent 
waves and inclusion of lossy materials. They also optimized parameters for PML. Others 
who did similar works were Cangellaris et al. [14], All above works focussed on split 
formulations. Later Chew and Weedon [18] introduced stretched co-ordinate approach, 
which provided .pathways to map PML in cylindrical and other co-ordinates systems. 
Also using stretched co-ordinate transformation, PML equations could be transformed to 
unsplit form, which is Maxwellian and more familiar. An alternative formulation for 
FEM, which was unsplit PML without co-ordinate transformation, was proposed by 
Sacks et al. [19]. Cangellaris et al. [17] and Gedney [20] adapted this unsplit formulation 
for FDTD. Gedney later also developed it for lossy and dispersive media [21]. He also 
analyzed patch antenna on isotropic substrate. Work of these researchers was for 
isotropic media. 

Garcia et al. [31] derived matching condition on a single interface between anisotropic 
medium and isotropic PML for an extra ordinary wave in 2-D space. However, it was of 
no practical use and could not be extended for 3-D space. Zhao et al. [32] proposed split 
formulation which they called Material Independent Perfectly Matched layer (MIPML) in 
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2-D dielectric space. They gave theoretical proof later [22] for this split formulation, 
which was very similar to what was proposed by Berenger, The proposed MIPML was 
used for 2-D arbitrary anisotropic magnetic media by Zhao [33] with almost similar 
approach by just replacing dielectric quantities with magnetic ones. Later Zhao et al. used 
the same MIPML for arbitrary dielectric and magnetic media for 2-D [23]. Zhao et al. 
[24] then extended this split formulation for 3-D but without any proof They used 
matching conditions similar to those of 2-D in three dimensions. 

In fact, it was Gedney’s work [20,21], which first time clearly hinted, and in his persona! 
communication to us, he confirmed that UPML could be used to absorb hybrid waves. 
Prompted by this, UPML has been adopted for our work on anisotropic media in this 
thesis. Using formulation for isotropic medium proposed by Gedney, It has been 
demonstrated that UPML equations are material independent if permittivity and 
permeability in Maxwell equations are combined with field quantities for replacing them 
with flux quantities. Then it has been shown numerically that this new UPML 
formulation matches Hybrid waves propagating in 3-D homogeneous and non- 
homogeneous arbitrary anisotropic mediums. For 2-D, a general proof has been 
developed for Split formulation to show that PML matches to anisotropic media. Zhao et 
al. [22] have also proposed a proof for matching PML to 2D anisotropic dielectric media 
using flux components. Proof in this thesis is without using flux components and it is for 
the media which is anisotropic in permittivity and permeability both. The proof has been 
done with field components only. In a publication Zhao et al. [24] have shown matching 
of 3-D anisotropic media to PML numerically but they have used split formulation which 
is quite memory intensive and secondly they have extended no theoretical proof In our 
work, it has been clearly shown that MIPML with unsplit formulation is far more 
efficient in terms of memory requirement. It has also been proved that extension of 
UPML to any media in 3D is possible, as UPML equations are inherently material 
independent. 

By adapting to MIPML, reflections from truncation boundary have been reduced. Its 
virtual reflection free performance gives highly accurate time domain results. By 
developing unsplit formulation, memory requirements have been reduced. By using 
Zhao’s advanced excitation scheme, source specification in FDTD code has been 
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simplified and reduction of errors in time domain results. This all together provides 
oscillation free frequency dependent characteristics of microstrip on an anisotropic 
substrate. 

Need for characteristics of planar circuits on arbitrary anisotropic substrate come from 
various reasons. One of the practical design difficulties of using isotropic substrates, such 
as alumina, is the significant variation in dielectric permittivity from different 
manufacturers or even from batch to batch from the same manufacturer. This essentially 
means that repeated measurements are required for accurate design of microstrip circuits. 
The use of anisotropic substrate with stable electrical properties, such as sapphire, 
alleviates this difficulty. Also sapphire has pot free structure and its flat surface makes it 
best known anisotropic substrate for MICs. But its anisotropic behavior requires new 
techniques to characterize microstrip devices on sapphire substrate. The TEM approach 
was used by Owens et al. [34] and Alexopoulos et al. [35] to determine Quasi-static 
characteristics of microstrip lines on anisotropic substrate with diagonal dielectric tensor. 
These characteristics are not accurate enough for high-speed pulse applications. In TEM 
approach dispersion is not considered which is significant at high frequencies. Mariki et 
al. [36] used TLM method to analyze enclosed microstrip on anisotropic substrate. Later 
Koike [37] et al. used Bergeron’s method in time domain. A unique problem with this 
method is that the dielectric interface and perfectly conducting strips are misaligned by 
half a space step. 

Our method is simple of all and takes dispersion into account. As stated earlier FDTD is 
efficient than other methods. MIPML is virtually reflection free. Advanced source 
excitation has been used so no spurious fields are induced. Use of new single cell method 
has removed whatever oscillatory behavior is there in the multiple cell method. These 
features make our results highly accurate compared to other methods. 

1.1 MIPML and Its Application to Present Work 

During recent years, great interests have been shown in using microstrip lines deposited 
on anisotropic substrates since substrate anisotropy could have important implications on 
the operation of microstrip circuits. With increasing complexity of geometry and material 
property, designing these circuits require more and more dedicated and sophisticated 
computer aided-design tools to predict the characteristics or performance of the circuits 
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as most of these staictures do not lead themselves to analytical solutions Therefore, the 
natural choice has been the FDTD method But in the past most of the simulation has 
focussed upon isotropic materials and also perfectly matched layers have been developed 
with independent derivations for TE and TM modes of propagation, assuming isotropic 
material. During the literature survey. It was discovered that all the deri\ations and 
theoretical proofs pertaining to perfectly matched layers have been done for TE and TM 
modes separately and then it has been assumed that in isotropic materials any propagating 
hybrid mode is superimposition of propagating TE and TM modes, so these derivations 
will hold good for propagating hybrid modes also. However, in a 3-D arbitrary 
anisotropic medium it is not possible to split any propagating hybrid mode into 
propagating TE and TM modes. So attempt had been made by us to carry forward these 
proofs to propagating hybrid modes within an anisotropic medium, as derivations and 
ultimate Maxwellian equations in PML for isotropic medium indicated that unsplit PML 
(UPML) should hold for anisotropic mediums also. Material independent nature of the 
coefficients of FDTD equations in UPML has been demonstrated, so it can be extended 
to anisotropic medium. By tiying the same matching conditions, which work for isotropic 
materials, numerically it has been prove that UPML works for anisotropic materials. This 
has opened possibility of using FDTD and PML combinations for circuits and radiators 
on anisotropic materials and predict their performance accurately. Here a generalized 
electromagnetic simulator ‘MIPML’ using FDTD and Unsplit Perfectly Matched Layers 
has been developed, which can analyze any material and geometry. This simulator is 
material independent and for each geometry and material, there is no need to rewrite 
whole algorithm. For each new geometry, one has to add boundary conditions for PEC in 
the working volume, specify dielectric permittivity tensor for all points in the working 
volume for the calculating electric fields from electric flux. In the algorithm of MIPML 
simulator, we deal with fluxes rather than fields. 

Towards the end of the thesis, dispersion characteristics of characteristic impedance and 
effective relative dielectric constant of microstrip have been worked out. Here both 
conventional multiple cell method [2] and new single cell method introduced by us, 
which is very effective in calculating above characteristics have been used. 
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1.2 Organization of the Thesis. 

Emphasis in chapter two and three will be on a brief introduction to FDTD method 
including a survey of it to demonstrate range of its applications and possible areas for 
future research. In fourth chapter, basic problems confronted by researchers in FDTD 
method have been taken up. In chapters, five, six and seven briefs details of Radiation 
boundary conditions and detailed discussion of material boundary conditions including 
split and unsplit approaches have been taken up. In eighth chapter, the scheme introduced 
in the present thesis, which has made material boundary conditions material independent 
is presented with its use in microwave circuits and antenna problems in anisotropic 
media. In nineth chapter, frequency dependent characteristics of microstrip on an 
arbitrary anisotropic substrate have been worked out. A new single cell method to 
calculate these characteristics has been introduced. Though the main objective of the 
thesis has been design of an electromagnetic simulator for microwave circuits and 
radiators in anisotropic media but in due course of evolution of algorithm it turned out 
that it encompasses material which are isotropic, and lossy. It can also cover dispersive 
and nonlinear materials with few very straightforward modifications. 
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Chapter 2 

FDTD Method: A Brief Survey 


Introduction 

There are three methods to predict electromagnetic effects in any geometry under 
consideration viz. Rigorous Analysis, Numerical Computation and Experimentation Of 
these methods. Numerical Computation is the newest and fastest growing approach Of 
the many approaches to electromagnetic computation, including MOM, FDTD, FEM, 
GTD, and Physical Optics, the FDTD technique is applicable to the widest range of the 
problems. As surveyed by K.L.Shlager and J.B. Schneider, there is sixty-fold increase in 
publications in 1996 as compared to 1985 and it is increasing day by day. The increase in 
publications is mainly due to various reasons. First can be attributed to increased 
computation facilities, which are growing cheaper very fast in terms of cost. Secondly as 
people are trying to encompass new problems of electromagnetics, they are recognizing 
new fundamental issues in FDTD and trying to address them. 

To select right approach for the problem at hand it is essential to study briefly about the 
existing research in that particular field. We have also done the same when we started 
thesis work. For the sake of completeness to the introduction to FDTD, we are providing 
a list of fundamental issues in FDTD. 

2.1 Fundamental Issues 

Following topics have been broadly recognized as fundamental issues of importance in 
the field of FDTD. 

2.1.1 Approaches 

There are two approaches to FDTD. 

(A) . Separate field formalism and 

(B) . Total field formalism. 

In first of these approaches the incident field is specified analytically through out the 
problem space and only the scattered field is determined. Here only scattered fields are 
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absorbed at the problem space outer boundaries. In the situations where incident fields 
are of much larger amplitude compared to scattered fields and only scattered fields are of 
importance, computing scattered fields only and absorbing them is advantageous from 
accuracy point of view. 

In the second approach total field is computed only over an interior subsection of the 
computational domain, while scattered fields are calculated in the remaining (exterior) 
portion of the grid. Radiation boundary conditions are applied over these scattered fields. 
To obtain this division of the computational domain into scattered-field and total-field 
regions, the incident field is specified over the boundary between these two regions. 

A comparison by Holland and Williams [3] demonstrated that, due to numerical 
dispersion, the total field is superior to scattered field approach. Secondly scattered field 
approach does not accommodate nonlinear media. However, for problems involving only 
linear medium and not containing shielded cavities, the scattered formulation is more 
desired approach. 

2.1.2 Absorbing Boundary Conditions 

These can be classified as 

(A) Differential-Equation based and other non material ABCs 

(B) Material ABCs. 

Differential-Equation based ABCs are approximation of one way wave equation [38] and 
Material ABCs are made up of artificial lossy material, which surrounds the computing 
domain. Both of these ABCs have been handled in details in subsequent chapters. 

2.13 Gridding 

There are various types of gridding schemes in use depending upon the problem under 
consideration. 

2. 1.3.1 Orthogonal Grids 

Cartesian gridding was the scheme proposed in original formulation but if any surface 
which varies smoothly, but not along the axes, then it has to be staircased. This 
approximation leads to significant errors in many problems. In addition, if any object 
under considerations has small-scale structure then in orthogonal gridding one has to use 
excessively fine gridding. 


8 



Therefore, it is natural that if any object is more naturally described in any other 
coordinate system other than Cartesian, it is always preferable to update equation for that 
particular system 

2.1.3.2 Subgridding 

Issue of fine structure is resolved by dividing problem space in sub-domains, which are 
gridded more finely than the rest of the problem. The key issue here is the coupling of the 
fine and coarse grids. Various interpolation techniques are used to define fields at the 
boundary of fine and coarse gridding. In addition to interpolation techniques, a 
discretised form of wave equation [4] is used to obtain fields on the boundary between 
the grids. 

2.1.3.3 Subceliular Techniques 

These techniques allow modeling of thin wires and similar objects without the use of 
excessively fine gridding [5] 

2. 1.3.4 Conformal Grids 

In this approach, the fields are expressed in terms of their covariant components (flow 
along a coordinate direction) and contravariant components (flow through a constant 
coordinate surface) and an integral formulation is used to obtain the update equations. In 
a nonorthogonal system, these covariant and contravariant are not collinear and auxiliary 
equations must be obtained to express one form in terms of the other. [6] 

2.1.4 Material Modeling 

Various types of materials with complex geometries are to be handled by FDTD. To 
solve these several techniques have been reported. 

2. 1.4.1 Frequency Dispersive Material 

One of the attractive features of the FDTD is that it is a wide band method and it will be 
incorrect if material parameters are assumed to be constant. To take care of them several 
schemes have been proposed. 
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2.1.4.2 Surface-Impedance Boundary Conditions 

These formulations have shown to provide significant computational savings over a full 
FDTD model. Several of SIBC formulations have been presented by researchers 

2.1. 4.3 Thin Material Sheets 

There are several schemes treating thin sheets of conducting and dielectric material but 
most of these treat only tangential components of the electric fields at air/dielectric 
boundary. A scheme suggested by Tirkas and Demarest [7] presented an FDTD model for 
thin dielectric sheets that treated not only the tangential components of electric fields at 
the air/sheet interface, but also the normal component of the field in the sheet. A 
comparison of the FDTD methods for modeling thin dielectric and conducting sheets was 
given by Maloney and smith [8] 

2. 1.4.4 Anisotropic Materials 

Treatment of these materials forms the core of this thesis. There have been several 
proposed schemes modifying original Yee’s algorithm for anisotropic materials but most 
of them focus upon the materials, which can be diagonalised. Recently, 1993 onwards, 
people have tried treatment of full tensor permeability and permittivity material in FDTD 
and our thesis tries to attempt extension of FDTD to these fUll tensor material and 
extension of PML to these problems. We have tried to use unsplit formulation of PML 
and though the treatment of these problems have been usually assumed to be quite 
memory intensive but we have tried a uniform algorithm through out the computing 
domain space and at the same time we have tried to conserve the memory requirement. 

2.1. 4.5 Nonlinear Materials 

Mostly in optics to model optical devices, techniques have been developed to treat 
nonlinear materials with FDTD. Apart from above, wherever material properties are 
function of fields present in media, several publications have been reported to treat the 
problem. Few to name are, nonlinear effects of ferrites, instantaneous nonlinearity of 
permittivity. Dispersive nonlinear Fabry-Perot cavities, propagation of famtosecond 
pulses in Dispersive and nondispersive nonlinear media, Lorentz dispersion, Kerr 
nonlinearities and incorporation of both Kerr and dispersive nonlinear effects to obtain 
temporal and spatial optical solitons from simulation of the full- wave time-domain 
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Maxwell s equations. The list may be endless but it is indication of prowess of FDTD 
method. 

2.1.5 Active and Passive Device Modeling 

This field has recently drawn attention and gained wide popularity. Tw’o approaches are 
used. In first analytical device models are used directly with the FDTD method. In 
second, lumped-element subgrid models are used in which device behavior is determined 
by other software. The second approach may be preferable in modeling of active devices 
with complicated equivalent circuit models. 

2.1.6 Transformations 

For complicated scatterers, time harmonic and pulsed illumination is used to determine 
the equivalent electric and magnetic currents on a virtual surface that completely encloses 
the scatterer. These currents are then transformed to the far field 

2.1.7 Digital Signal-Processing Techniques 

By using digital signal processing techniques in conjunction with the FDTD method 
relevant data such as frequency-domain scattering parameters can be extracted from 
FDTD simulations of the shorter duration than would otherwise be possible. 

2.1.8 Techniques to Reduce Numerical Dispersion Errors 

To reduce numerical dispersion errors many people suggested higher order 
approximations of space and time derivatives. But ,recently proposed techniques ,such as 
those employing wavelets, make a significant departure from traditional methods in order 
to reduce dispersion errors. One such technique, derived from theory of wavelets, is 
Multiresolution time domain (MRTD)[9]. To suppress spurious modes, researchers have 
used Daubechies wavelets. 

2.1.9 Radiating Structures 

Initially FDTD was used to model scattering from objects. However by including the 
feeding source within the grid , it is possible to model radiating structures. This has led to 
the use of FDTD to model antennas. Papers have been published on following 
(A) Simple antennas like Cylindrical and- conical monopoles, thinwire dipoles and thin 
wire Yagi antennas etc 
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(B) {lorn Antennas 

(C) Antennas for Pulse Radiation 

(D) Microstrip antennas 

(E) Dielectric Resonator Antennas 

(F) Hand held Antennas 

(G) Antenna Arrays 

(H) Other Radiating surfaces like cavity backed antennas. Active antennas etc. 

2.1.10 Microwave Devices and Guiding Structures 

The FDTD has found wide applications in modeling microwave de\ ices and guiding 
structures. In this area, several techniques to improve the efficiency of the FDTD have 
been developed. Such techniques include the compact two dimensional FDTD method for 
guided-wave structures, ABCs specifically developed for Guided wave problems, and the 
previously Discussed advanced schemes for gridding and reducing numerical dispersion. 
To name a few areas, we are giving a list 

(A) Waveguides, Feeds, Junctions, and Resonators 

(B) Microstrip 

(C) Vias , interconnects and transmission lines 

(D) Algorithm improvements 

There has been active research in reducing the computations needed to analyze three- 
dimensional guided-wave structures. These approaches utilize compact two-dimensional 
FDTD method, the MRTD method, and conformal FDTD methods. The compact two- 
dimensional FDTD method is based on introducing a phase shift along the direction of 
propagation, which reduces three dimensions to two dimensions 

2.1.11 Discrete Scatterers 

Research has mostly been confined to calculation of Radar cross section (RCS of 
scatterer). In addition to that, modeling of perfectly conducting or dielectric scatterers has 
been done. 
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2.1.12 Infinite and Periodic Structures 

These include periodic guided wave structures, periodic phased array antenna, metallic 
strip-grating etc. 

2.1.13 Ground Penetrating Radar 

Prompt by the use of FDTD to study the scattering from a discrete object or the radiation 
from a source in the presence of a planar boundary, researchers have tried problems 
related to subsurface radars and ground penetrating radars. 

2.1.14 Hybrid Techniques 

The FDTD method has also been combined with other techniques where it would be 
difficult or infeasible to model the complete geometry in a single space lattice. A hybrid 
FDTD/method of moments approach [10] based upon schelkunofTs field-equivalence 
theorems to investigate coupling problems involving loaded cavities have been reported 
long back. Others are hybrid spectral/FDTD, Hybrid ray/FDTD, Hybrid FDTD/partial- 
eigen function-expansion method and hybrid FDTD asymptotic method etc. 

A Web site http ://www. eecs. wsu . edu/ ~schneidi/ fdtd-bib. html presents listing of papers 
on FDTD. Here list is exhaustive and very useful for finding references on any FDTD 
related topics in which someone wishes to start his research career. 
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Chapter3 


The FDTD Method 

Introduction 

In this chapter, we are going to review the basic Maxwell's equations and FDTD 
algorithm based on these equations We have focussed upon central differencing scheme, 
using the same to discretise Maxwell equations and construction of FDTD algorithm for 
full wave analysis of a most generalized problem. Also we have tried to include 
algorithm suitable from programming point of view. 

3.1 Maxwell Equations 

The Maxwell curl equations in Cartesian Coordinates are written as 
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Where 





Ex, Ey, Ez, Hx, Hy and Hz are the Cartesian components of electric and magnetic fields, 
s is electrical permittivity and p is magnetic permeability, 
p is magnetic resistivity and a is electrical conductivity. 
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3,2 The Central Difference Scheme 

Using following scheme of discretisation, above equations can be solved Central 
Differencing scheme has been used because it is second order accurate 
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3.3 FDTD Algorithm for Full Wave Analysis 

If, we discretise Equation (3.4) as per above scheme. The resulting equation is 
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As evident from above equation, that dealing with half indices is problematic, as most of 
the programming languages do permit matrix indices to be integers only. If Yee’s [1] 
original cell is modified as shown below then above problem can be easily circumvented. 
Here position of field components is kept same, only the way they are addressed is 
different. Here indices are associated with each cell and all fields components in this cell 
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difterent. Here indices are associated with each cell and all fields components in this cell 
have the same indices though their position is implied in the cell and is different from 
each other. This means in a cell, field components are still half space step apart from each 
other but their indices are same as shown below 



F igure3 . 1 YEE' s cell 


As is evident that above formulation makes programming an easy affair with no 
botheration of half indices. Above scheme of Yee is called ‘leapfrog in time’ scheme. 
Here first H is evaluated at time (N+1/2) from H at (N-1/2) and E at N. Then E 
components are evaluated at time (N+1) using H at (N+1/2) and E at N. Again N is 
increased By 1 and same procedure is repeated. This interleaves E and H temporally and 
results in a centered difference or ‘ leapfrog in time’ approach. Following set of equations 
result by using above scheme. We are giving a complete set of equations for all six 
components 

H Magnetic field, p Magnetic resistivity 
£ Electric field, a Electrical conductivity 
s Permittivi ty, p permeabili ty 
Ar,A>’,Az space steps. A/ time steps 
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( 3 . 18 ) 


3.4 FDTD Code Architecture 

At this Stage it is apt to consider FDTD code requirements and architecture of the 
program. Inefficiently written program may take large memory and execution time may 
become so large that testing and debugging may not be possible. Here we will consider 
FDTD code architecture for Perfectly Matched Layers .It may slightly differ from 
application to application but in general, it is applicable to almost all applications. 

The code architecture is summarized as; 

• Driver 

• Problem space setup 

• Test object definition 

• E, H field algorithms 

• Absorbing Boundary Conditions 

• Data saver 

• Data processing for the given problem 

• Results display 

3.4.1 Driver 

• Calls problem setup subroutine and the test object definition subroutine 

• Time steps over index N 
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• While looping over N, calls h. H subroutines and Absorbing Boundar>’ Condition 
subroutines 

• At appropriate time steps or at each step it calls data saver to store sampled data 

• After completion of looping over index N routines for data processing are called, 
results are computed and displayed. 

3.4.2 Problem Space setup 

• Sets the problem space size 

Sets the number of cells in each dimension 
Sets the cell size (Ax, Ay, Az) 

• Calculates the At time step according to courant stability condition 

• Calculates the constant multipliers 

3.4.3. Test Object Definition 

Cells or individual field components in the cells are flagged indicating their composition. 
This dictates how the E, H algorithm is to process the data: as perfect conductor, lossy 
dielectric, free space, or other more complicated material. It is usually convenient to set 
default material to free space; the array of flags can be read for a preprocessing check of 
the geometry and composition of the object. 

3.4.4 E, H Field Algorithm 

Calculate the response of a component from its own prior time value and that of the 
nearest-neighbor field quantities (Es around Hs and Hs around the Es) according to the 
type of material present at that component location: 

Free space 
Lossy dielectric 
Lossy magnetic 
Perfect conductor 

3.4.5 Absorbing Boundary Conditions(ABCs) 

These conditions truncate computing domain, ideally without causing any reflection of 
fields from the truncation boundar>'. Again depending upon ABCs, implementation may 
, differ. ABCs can be split-up in several regions as face ABCs, Edge ABCs and Comer 
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ABCs [ 1 7]. For each region there will be a different algorithm and they will be computed 
one after another. In ABCs, equation constants may he different from working volume 
and the>' may be different in different ABC regions. So if they are treated one after 
another, code may be quite long though it maj' take less memory and computationally it 
may be fast. But obvious disadvantages are that it is quite time consuming to develop this 
kind of code and once code is completed it is a big headache to debug it due to the 
enormous size of code and a number of subroutines. In our opinion if sufficient memory 
is available than treat whole computing domain as made up of single ABC region and 
define value of elements of a constant multiplier matrices, depending upon the region 
they represent. This makes memory requirement large but code is uniform, very compact, 
easy to debug and execution time does not increase significantly. In fact there is no need 
for separate ABC algorithm. The E, H field algorithm takes care of ABCs. 

3.4.6 Data Saver 

Saves response data such as E and H field components, currents, or other quantities in the 
FDTD computation space in arrays at chosen time steps. 

3.4.7 Data Processing and Result Display 

This basically consists of calling library routines for Fourier transforms, various 
mathematical operations and graphic displays of the fields and writing results to a file for 
later reference. 

A generalized flow chart for any FDTD code is given in figure (3.2). 
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FDTD FLOW CHART 

Figure (3.2). 
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Chapter 4 

FDTD Fundamentals 


Introduetion 

In this Chapter, the practical difficulties for implementation of FDTD method have been 
addressed. There are three main issues viz. Cell size. Time step size, and Source 
specification. 

4.1 Cell Size 

There are two basic factors, which affect cell size. As per the Nyquist sampling theorem, 
highest frequency in the signal must be sampled at least twice. So X = 2 Ax gives the 
maximum cell size which could be taken. For accurate results usually a smaller fraction 
of wavelength is taken .The limit varies depending upon the accuracy demand and 
usually it is XyiO to X/20 or even smaller. Such a high level of sampling is required 
because at any particular time step the FDTD grid is a discrete spatial sample of the field 
distribution. From the Nyquist sampling theorem, there must be at least two samples per 
spatial period (wavelength) in order for the spatial information to be adequately sampled. 
Because our sampling is not exact, and our smallest wavelength is not precisely known, 
more than two samples per wavelength are required. Another related consideration is grid 
dispersion error. Due to the approximation inherent in FDTD, waves of different 
frequencies will propagate at slightly different speeds through the grid. This difference in 
propagation speed also depends on the direction of propagation relative to the grid. For 
accurate and stable results, the grid dispersion error must be reduced to an acceptable 
level, which can be readily accomplished by reducing the cell size. By looking at 
Numerical Dispersion relation [40], it is quite evident that as spatial step reduces 
Numerical Dispersion relation reduces to Analytical Dispersion relation. 

Another cell size consideration is that the important characteristics of the problem 
geometry must be accurately modeled. Normally if cell size is A/ 10 to XIIQ or even 
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smaller then this requirement is usually met but in certain cases if it is not so, then finer 
gridding is required. This sometimes leads to a situation where grid size is excessively 
large and some times object can not be properly modeled by finer grid e.g. thin wires. In 
such cases Sub-cellular techniques [5] are used Many a times subgridding [4] comes 
handy. Depending upon nature of the problem gridding is done and size of computational 
domain is determined. 

4.2 Time Step Size 

Once the cell size is determined, The maximum size of the time step At immediately 
follows from the Courant condition [41]. The basis for the Courant condition is that in 
one time step any point on a plane wave must not pass through more than one cell. It is 
because during one time step FDTD can propagate the wave only from one cell to its 
nearest neighbors. In the direction perpendicular to lattice planes, wave propagates most 
rapidly. So with this we can say 

V* At < Ax for stability in one dimension (4.1) 


For three dimensional space 


V Ai< 



(4.2) 


When equality holds, the discretised wave most closely follows the actual wave 
propagation and grid dispersion errors are minimized. Smaller values of At does not give 
most accurate results. However exception to this occurs in case of conductive materials 
and nonlinear materials where time step size is smaller than what is required for equality. 

4.3 Source Specification 

The type of source to be used is dependent upon the problem under specification but 
there are several popular types of sources, which have emerged as experience has grown 
with FDTD method. 
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4.3.1 Compact Pulse 

SourceiJ ) = (l 0 - 1 5 cos + 6cos4;^/-cos bitl ) for t < 1 ns 

Sourceit) = 0 for t > 1 ns (4.3) 

Use of compact pulse is reserved for the applications where response over a wide band is 
not required and grid dispersion errors are to be kept to a minimum. Examples are testing 
or comparison of various Radiation Boundary Conditions. This compact pulse has very 
little high frequency contents and thus dispersion errors are hold to a minimum. 

4.3.2 Gaussian Pulse 

Here we take ‘pw’ as half width of the Gaussian pulse in seconds and the width is defined 
as the distance between the points where pulse reduces to 5 % of its maximum value, ‘x’ 
is the location of the peak of the pulse in sec. 


r = pw 


alpha = 


Sonrce{t) = e 


0. 2S* pw ) 

^ Lalpha*{t-rf) 


(4.4) 


Maximum frequency content of Gaussian pulse is approximately given by 

/max “ (4-5) 

2* pw 

This source is preferred for the applications where response over a ■wide band of 
frequencies is desired. Its Fourier transform is also Gaussian. Whenever space is filled up 
with penetrable material then spatial cell size is calculated keeping in mind the 
wavelength of the wave in that material. So accordingly cell size reduces by a factor of 
square root of dielectric constant which leads to reduction in value of time step almost by 
the same factor and if pulse width is specified in terms of time steps its width reduces to 
half in terms of time and frequency content goes double which in turn will lead to a lot of 
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grid dispersion error. So care should always be taken to specify pulse width always in 
seconds and this reduces chance for making any mistake in source excitation. 

Another guideline for using the Gaussian pulse is to take its Fourier transform and look 
for a 20 dB down frequency. If at this frequency spatial step size is one tenth of the 
wavelength than we can have accurate results upto that frequency. As far as noise and 
stability are concerned, if we have appreciable energy in high frequency components then 
it could be a problem. But by taking tau sufficiently large this possibility of high 
frequency components in pulse is reduced. Also to put quantitatively, we have put 
another criterion that at wavelengths four times the minimum cell size, pulse amplitude 
should be down by 120 dB. 

4.3.3 Sinusoidal Wave 

Wherever response for a single frequency is required, sinusoidal source is used. 

Source{t) = sin(2 * tt* f *t) (4.6) 

One major issue in source consideration is the point of excitation and mode of excitation. 
Usually excitation in dominant mode is preferred but some times it is not known 
accurately. Reason for exciting in dominant mode is that, less number of other unwanted 
modes are generated. If it is not known accurately then closest approximation is done, so 
less number of other unwanted modes are generated and little distance is required for 
them to travel to die down and what propagates ultimately down the line is the dominant 
mode. Dominant mode leads to reduction in size of computing domain. 

In addition, point of excitation is very' important. In paper by Zhang et al. [2] the 
excitation was done on front wall of a computing domain inside of which, microstrip 
discontinuities were analyzed. They excited microstrip on the front wall by setting 
tangential electric field under the strip as gaussian pulse and rests of the points on plane 
were set to zero. This amounted to an electric wall condition, which induced a dc current 
or tangential magnetic field on the front surface and near by. When reflected pulse 
arrived they switched on an absorbing boundary condition. If boundary conditions are 
applied on this local magnetic field, errors accumulate and solution blows up. So after 
pulse left the plane and before it reflected back from discontinuity, boundary conditions 


26 



steps in side the computing domain to avoid inclusion of this local DC field. In summary 
this is a very complicated arrangement and size of computing domain also increases 
Another scheme used by sheen et al. [11] was to make source plane a magnetic wall with 
tangential electric field existing on it. Magnetic wall condition was imposed by them 
using image theory. They set Ht.™ outside the magnetic wall equal to -Htan inside the wall 
in terms of amplitude and, then remaining electric field components on the source wall 
were calculated using the finite-difference equations. This gave little distortion of the 
source. 

A very novel and highly accurate scheme was proposed by Zhao et al. [27,28]. In this 
scheme source plane was put few cells inside the computing domain at near-end terminal 
as shown in figure (4.1) and plane wave source condition [41] was applied. Though the 
excitation was done beneath the microstrip but remaining points on the plane were 
calculated using normal FDTD algorithm and at excitation points following equation was 
applied to simulate plane wave source condition 



Source plane 


Figure 4.1 Advanced source excitation scheme 
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1 .Above equation is for the points lying below the strip in the plane inp. For rest of the 
points in the inp plane simple FDTD equations are applied 

2. Above equation assumes that all conductivities and resistivities are zero. For lossy 
materials the equation need to be modified accordingly. 

Within the region 1, the EM fields, contain the incident wave (propagates in +y direction) 
and reflected wave from any discontinuity, if present. Where as the EM fields in region2, 
contain the incident wave (propagation in -Y direction) and the reflected wave. In region 

2 the incident wave propagating in -Y direction and reflected waves are immediately 
absorbed on near end terminal plane, while the incident wave propagating in +Y direction 
is used to examine the discontinuous system itself 

It is also to be noted that while exciting in three dimension space, using a point source, 
above equation is not to be used. It has been observed that wave does not propagate. Only 
last term of incident source is to be used and rest of the terms are dropped in this case. 
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Chapter 5 


Radiation Boundary Conditions 


Introduction 

There are two broad categories of the problems encountered in EM theory In one type of 
problems. Initial and Boundary values are known whereas another category includes the 
problems where boundary conditions are not known. These problems are called open 
problems. When these problems are dealt by FDTD method the natural question is ‘how to 
define computing domain’. As whole infinite space can not be taken as computing domain 
because of limitation of computer memory, we need to keep it limited. So computing domain 
is limited to a size where it can accommodate structure or geometry to be analyzed. Once 
computing domain is limited then obvious problem arises that at the boundary of this domain, 
how to apply FDTD equations, as central differencing schemes need, value of fields inside as 
well as outside of the boundary. Since we do not have values outside the boundary, we can 
not use simple FDTD equations here. Some new set of equations is required which simulate 
propagation of waves to infinity. These equations are expected to absorb whatever waves 
strike at boundary without giving any reflection. Since properties of waves arriving at the 
boundaries is not known before hand, designing of these absorbing equations, which are also 
called absorbing or radiation boundary conditions (ABCs or RBCs), is a difficult task. 

Primary requirement of an ABC is that irrespective polarization, angle of incidence and 
frequency of incident waves, it must absorb all of them without any reflection. This reflection 
free boundary will permit the FDTD solution of the fields inside the computing domain to be 
true at all time steps. Though the ABCs described here have been rendered obsolete because 
of the invention of Perfectly Matched Layers, which provides performance that is much more 
superior. This material is still important from the point of view of academic interest and for 
implementing these ABCs where computer resources are extremely limited. 

In this section, the theory and numerical implementation of one way wave equation ABCs, a 
very useful class of radiation conditions in Cartesian coordinates has been presented. One of 
them, the Mur’s second order radiation boundary condition, is appropriate for effectively 
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truncating a two- or three- dimensional FDTD space lattice with an overall level of spurious 
reflections of l%-5% for outer lattice planes located 10 - 20 space cells from a target surface 

5.1 One Way Wave Equations 

A partial differential equation, which permits wave propagation only in one direction, is 
called a “one way wave equation.” Figure (5 1} below shows a finite, two-dimensional 
Cartesian domain, on which the time dependent wave equation is to be simulated. In the 
interior ofO, a numerical scheme which models wave propagation in all direction is applied 
On rQ, the outer Boundary of Q, only numerical wave motion that is outward from Q is 
permitted. The boundary must permit outward propagating numerical wave analogs to exit Q 
just as if the simulation were performed on a computational domain of infinite extent A 
scheme, which enacts a one way wave equation on dQ for this purpose is called a Radiation 
Boundary Condition. 

Y. 

dQ 

— X 

Fig 5.1 Numencal plane wave analog incident upon left grid boundary' of 
a 2 d Cartesian comDutational Domain 



5.2 Derivation By Wave Equation Factoring 

The Derivation of a RBC whose purpose is to absorb numerical waves incident upon the 
outer boundary of a finite difference grid can be explained in terms of operator factoring. 
Consider a two dimensional wave equation in Cartesian coordinates 

( 5 . 1 ) 

c 

This is a compact representation of following wave equation 

d-U d-U 1 d'U 
ex' dy' c' dr 
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Where U is a scalar field component; the subscript xx, yy, and tt denote second partial 
derivatives with respect to x, y, and t, respectively; c is the wave phase velocity The partial 
diflerential operator here is 


L = D;+d;~-Ld; (S 3 ) 

C'' 


Which uses the notation 


CdC 
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a-= 


The wave equation is then compactly written as 

LU = 0 (5.5) 

The wave operator, L, can be factored in the following manner: 

LU = L"L-f/ = 0 (5.6) 


Where 


L (5.7) 

c 


With 



(5.8) 


The operator L*, is similarly defined except for a “+”sign before the radical. 

Enquist and Mazda [ 38 ] showed that at a grid boundary, say at x= 0 , the application of L‘ to 
the wave function, U, will exactly absorb a plane wave propagating toward the boundary at 
an arbitrary angle, 0 . Thus, 


L-U = 0 (5.9) 



If applied at x=0, functions as an exact analytical RBC which absorbs wave motion from the 
interior of the spatial domain, Q. The operator, L*, performs the same function for a plane 
wave propagating at an arbitrary angle towards the other x boundary in Figure (5. 1 ) at x = h. 
The presence of the radical in equation (5.7) classifies V as a pseudo-differential operator 
that is non-local in both the space and time variable. This is an undesirable characteristic in 
that it prohibits the direct numerical implementation of (5.9) as a RBC. .Approximation of the 
radical in (5.7) produces RBC’s that can be implemented numerically and are useful in FDTD 
simulations with minimum reflections over a range of incident angles The RBC, used in 
FDTD electromagnetic wave codes, is simply a two term Taylor series approximation to the 
radical in (5.7), given by 

(5.10) 


Substituting (5.10) in (5.7), we obtain 
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Multiplying (5.11) through by Dt, and identifying the differential operators as partial 
derivatives, we obtain the following analytical RBC which can be numerically implemented 
at the X = 0 grid boundary 

(5,12) 

c 2 ■■ 


Equation (5.12) is a very good approximation to exact RBC of (5.9) for relatively small 
values of S, which satisfy the Taylor series approximation of (5.10). This is equivalent to 
saying that (5.12) presents a nearly reflection less grid truncation for numerical plane wave 
modes which strike the x=0 grid boundary at small values of the incident angle, 0. Analogous 
approximate, analytical RBC’s can be derived for the other grid boundaries 

(5,13) 

C 2 " 
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at X = h 




U,; +-(!„ (5.15) aty = h 

c 2 

For the FD-TD simulation of the vector Maxwell's equations, the RBC’s of (5. 12)-(5. 1 5) are 
applied to individual Cartesian components of E or H that are located at, and tangential to. the 
grid boundaries. 

The Derivation of RBC’s for the three-dimensional case follows the above development 
closely. The wave equation, given by 

u^+u„+u^-\u.^0 (5.16) 

c 

has the associated partial differential operator 

L = Dl+D;+D;-^D; (5.1?) 

c 

L can be factored in a manner of (5.6) to provide an exact radiation boundary, L" having the 
same form as that of (5.7), but with S given by 



Again, L”, applied to the scalar wave function, U, at the grid boundary will exactly absorb a 
plane wave propagating toward the boundary at an arbitrary angle. 

Using the Taylor series approximation of (5.7), we obtain an approximate RBC at x = 0 in 
differential-operator form 



£l 

c 



2D, 



7 = 0 


(5.19) 
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Multiplying (5.19) through by Dt, and identifying the differential operators as partial 
derivatives, we obtain the corresponding approximate, analytical RBC which can be 
numerically implemented at the x=0 lattice boundary 

(5.20) 

c 2 "2 

Equation (5.20) is a very good approximation of the exact RBC of (5.9) for relatively small 
values of S given by (5.18). This is equivalent to saying that (5.20) presents a nearly 
reflection less lattice truncation for numerical plane wave modes which strike the x = 0 lattice 
boundary close to broadside. Analogous approximate, analytical RBC’s can be derived for 
the other lattice boundaries 

2 f 

(5.22) aty = 0 

c 2 2 

f'v, (5.23) at y = h 

c 2 2 

- -U„ = 0 (5.24) at z = 0 

c 2 2 '* 

(5.25) atz=h 

For the FD-TD simulation of the vector Maxwell’s equations, the RBC’s of (5.20) - (5.25) 
are applied to individual Cartesian components of E or H that are located at, and tangential to, 
the grid boundaries. 

Equations (5.12)-(5.15) for two dimensional grid, and equations (5.20)-(5.25), representing 
three-dimensional grid have been found to be very effective if differencing scheme suggested 
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by Mur [39] is used. Mur’s Scheme has been shown lo truncate lattice space with just 1% - 
5% reflections for arbitrary targets. 

5.3 Mur Differencing scheme 

Consider the following Figure (5.2) used for illustrating Mur’s RBC 


WO.ifl) 


Y=fj+| |,\\ 



Y- j \\ 




WJ-l) 

Y=(j-l)Ay 

X=0 

X=-iX 



Fig5.2. Points near the xK) Boimdaiy used in Mir diflfererrar^ 
Sdieme 


A Finite difference based scheme for 2-D RBC equations (5.12)-(5.15) and 3-D RBC 
equations (5.20)-(5.25) was introduced by Mur [39], We are illustrating the scheme for 2-D 
grid case at x = 0. Mur applied equation (5.12) at point W (l/2,j) and used central 
Differencing to write partial derivatives. The mixed partial x and t derivatives on left-hand 
side of equation (5.12) are written out using central differences 


cW 






cbc 




dK 


1 

2 y 


.;.n 


2A/ 


’w"^’(i,;)-w"^'(o,y)' 


"W" ‘(l,i)-W"-'(0,y)' 

Ax 


Ax 


2.A/ 


(5.26) 
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Next, partial t derivatives on the left hand side of (5.12) is written out as an average of time 
derivatives at the adjacent points (0,j} and ( l,j) 






rf cl J 

2 [ A/ 

. ^^''"'(l 7)-2»-"(l,7) + lf'’ '(1. 7)1 

Ar ~ ’ i 


(5 27) 


And, the partial y derivative on the left hand side of (5.12) is written out as an average of y 
derivatives at the adjacent points (OJ) and (l,j) 




aV". .X d'w” . . 

-TT-(0’^)+^rr-(i,7) 

cy ^ 

W''{0J + \)-2W''{0j)+W”{0J-l) 
A/- 

W”{\J + l)-2W’'{lj)+W"il,j-iy 
Ay~ 


+ 


(5.28) 


Substituting the finite difference expressions of (5.26)-(5.28) and solving for W""’(0,j),we 
obtain the following time-stepping algorithm for components of W along the x = 0 grid 
boundary which implements the Taylor’s series of RBC of (5.12) 

cAt + Ax 

+ - 7 ^[(f(i,;)+(f"(o,j)l 

cA/ + Ar 

+ +l)-2(f' (0,y)+H'"(0,j - 1) 

+ H'-(l,; + l)-2(f(l,y)+(»'-(l,./-l)l (.‘i.29) 

Analogous finite-difference expressions for the Mur RBC at each of the other grid 
boundaries x = h, y=0, and y = h, can be derived by proper substitutions into (5.13), (5.14), 
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and (5 15), respectively, in the same manner. More simply, these Mur RBC can be obtained 
by inspection from (5.26)-(5.29) using coordinate symmetry arguments. 

RBC for 3-D case can be derived, if Figure (5.2) is now supposed to be representing 
individual Cartesian components of E and H located in lattice plane z - kAz at x = 0 lattice 
boundary. Here, the Mur scheme involves implementing the partial derivatives of (5.20) as 
numerical central differences expanded about the auxiliary W component, Wn (l/2,j, k), 
located one half space cell from the grid boundary at (0,j, k). The partial derivatives. W^i, Wn. 
and are identical in form to (5.26), (5.27), and (5.28), respectively, and are evaluated in 
lattice plane z = kAz. The partial derivative, Wzz, is expressed as an average of z derivatives 
at the adjacent points (0, j, k) and ( 1 , j, k) 




~ 2 


+ 


C'Z' cz 


W'-{0,J,k + l)-2W''{0J,k)+W'’{0J,k-\) 
W’'{\J.k + \)-2W''{\J,k)+W''{\J,k-\) 


(5.30) 


Substituting these finite-difference expressions into (5.20) and solving for W"’^’(0, j, k), we 
obtain the following time stepping algorithm for components of W along the x=0 lattice 
boundary which implements the Taylor series RBC of (5.20). 




cA/ + Ar 

cAt + Ax 
2Ay{cAi + Ax) 

+ IV{l,J + lk)-2W(}Jj)+IV(lJ-i.k)] 

+ . [^•{0j^t + l)-2W-{0J.k)+IV-(0J.k-\) 

2Az~{cAl +Ax) 

+ W{\J,k + ])-2W''{\J,k)+W'‘{lJ,k-\)] ( 5 . 31 ) 

Analogous finite difference expressions for the Mur RBC at each of the other lattice 
boundaries, x = h, y = 0, y = h, z = 0, and z = h, can be derived by substituting into (5.21)- 
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(5,25) Respectively, in the same manner More simply, these Mur RBC’s can be obtained b\ 
inspection from (5.26)-(5 28) and (5 30) using coordinate symmetry arguments 

5.4 Special Corner RBC 

Upon inspecting (5,29) and (5.31), it is clear that the Mur’s finite difference scheme for the 
two terms Taylor series RBC's cannot be implemented for the field components located at 
grid corners. It is due to the fact that some of the necessar\' field data used in the Mur 
expressions at these points are outside of the grid and not available It is necessary to 
implement a special corner radiation boundary condition at these points which (1) utilizes 
available field data in the grid; (2) yields acceptably low levels of reflection of outgoing 
numerical wave modes; and (3) numerically stable. 

Figure (5.3) illustrates the two-dimensional grid geometry for a simple and stable special 
comer RBC used successfully for a variety of two and three-dimensional FDTD simulations. 



Figure 5.3: Points near the x=0,y=0 grid comer used in the special 
comer radiation boundary condition (square grid case). 


Here basic assumption is that wave propagates along a radial line from the center of a 
domain. So field at comer point W (0,0), is taken to be just the time retarded value of an 
interior field, W, located along a radial line connecting the comer point to the center of the 
grid. Also it is assumed that cAt = dll is maintained, so that if, W is located exactly one cell 
width, 6, inward along the radial line. The time retardation of the outgoing numerical wave in 
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propagating from W to VV (0,0) is exactly two time steps. Over ail, the special corner RBC is 
given by 


IF'"' 


-.L.ujy 


n- 


(5.32) 


Where t raiii.ii is the attenuation factor for the radially outgoing wave .In two dimensions from 
figure (5.3) 


f. 


radial 



(5.33) 


W" ' - (l -sinaX^ "COsa)PF" ‘(0,0) 

+ (l - sin a)cosa W ' (l,o) 

+ sin C!r(l - cosor)fF”' (0,l) 

+ sin a cosorPF"''(^--0 (5-34) 


Where dceirter is the radial distance, in cell-widths, from W to the center of the grid, and a is 
the azimuthal angle of the radial line at W (0,0). Note that the value of W is determined by 
simple linear interpolation of the four surrounding field values including W (0,0) at time step 
n-1 . Extension to three dimension is straight forward, yielding for W (0,0,k) 


f, 


radial 




(5.35) 


fV” ' = (l - sin - cos sin aX^ ” P cosotX^" ' A^) 

+ (l - sin py^ - cos P sin a)cos P cosalV " ' (l,0,A^) 

+ (l - sin y9)cos p sin a(l - cos P cosa)fV (0,1,^) 

+ (i -sin yff) COS' sin a cosa.fr (l,I,^) 

+ sin p(l - cos P sin aXl - cos P cosa)fr (O.O, /: + 1) 
+ sin y3(l -cos>9sina)cos/?cosafr (l,0,^ + l) 

+ sin p cos P sin a(l - cos P cosa)ff (0,1, A + 1) 

+ sin Pcos~ psinacosaW" '(l,l,^ + 1) 


(5.36) 



p is the azimuthal angle of the radial line at W fO.O.k)Note that the value of W ’ is 
determined by simple linear interpolation of the eight surrounding field values, including 
W(0,0,k), at time step n-1. Other corners can also be worked out in a similar manner or 
simply by inspection equations can be written using coordinate symmetry argument 
Another corner RBC, which also works well, is given as below 

W"*Uo,0) = W"m) + - M''"{0,0)) (5.37) 

V * deli -i-del 


We have worked with both corner RBC and results are not very different 

5.5 Generalized and Higher-Order RBC’s 

Trefethen and Halpem [42] proposed a generalization of the two terms Taylor series 
approximation to the radical in equation (5.7), considering the use of the rational function 
approximation 



(5,38) 


On the interval [-1,1], where pm and qn are the polynomials in S of degree m and n, 
respectively; and r (S) is said to be of type (m, n). With S = c Dy / Dt, the [-1,1] 
approximation interval on S is equivalent to approximation of the exact one way wave 
equation of equation (5.9) along the x = 0 grid boundary for the range of incident wave angles 
0 = -90° to 0 = 90° 


For example, by specifying r (S) as a general (2,0) approximant, the radical is approximated 
by an interpolating polynomial of the form 





(5.39) 


Resulting in the general second-order, approximate, analytical RBC, 

U.-—U.-p.cU„=0 (5.40) 

C 

The choice of the coefficients, po and p 2 , is determined by the method of interpolation that is 
used. Standard techniques such as Fade', least-square, or Chebyshev approximation are 
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applied with the goal of interpolating the radical optimally over the 1-1,1 ] range of S, there by 
producing an approximate RBC whose performance is good over a wide range of incident 
wave angles. Mur’s two term Taylor series approximation of (5.9) is now seen in a more 
general sense as a pade’ (2,0) interpolation i.e. with coefficient s po ■= 1 and p2 = -1/2 in 
(5.40). Higher order approximation to the radical in equation (5 7) gives better ABCs, for 
example, the use of the general type (2,2) rational function 

Gives the general third-order, approximate, analytical RBC 

qJU„ + - p,cV„ = 0 (5.42) 

c 

Appropriate selection of the p and q coefficients in (5.42) produces various families of RBC’s 
as suggested in [42]. The only problem is that as higher order approximations are used, 
discretized equations become quite large. Third order equations may contain as many as 50 
terms or more at each of the grid truncation boundary. These also require increased storage. 
Since the invention of Perfectly Matched Layers, this third and higher order RBC’s have 
become obsolete. Only Mur’s RBC is used because of its simplicity and that too, in situation 
where demand on accuracy is not very high. 

5.6 Numerical Experimentation with MUR’s RBC 

We had planned to use Mur’s RBC when we started with the thesis work as till that time we 
were not aware of PML or its MIPML variant. That time we had developed code SO ABC for 
second order Mur RBC for use later in project. Our results were in excellent confirmation 
with results produced by Berenger in his paper [13], We have performed test used in [13], 
We have considered the same 100*50 cell domain, with the same space and time increments, 
respectively 1.5 cm and 25 picosec. The domain was surrounded by MUR RBC. Our 
computation was done for both TE and TM case .In TE case excitation was done at magnetic 
point and for TM case excitation was done at electric point. Excitation was done using a 
compact pulse 
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(5,43) 


H , ( 50 ^ 20 )^ ~(10-15cos2;z/4 6cos 4;7/ - cos6;z/) 

H ^{50^20)^0 ifi>lns 

A reference solution was computed using a large domain of 400*400, allowing a boundary 
free solution during a 500 time steps. Denoting Hz (ij) the field in the test domain and 
(i,j),its counterpart in reference domain ,two kinds of results are presented and compared with 
those of Berenger[l 3], 

First the reflected field along a boundary (j = 1), at the time step 100, normalized to the peak 
value of the field at point (50, 1) 

R(i) = (5.44) 

And, second, the 1? norm of the error on the 100*50 domain as a function of time, 

too 50 - 

(5-45) 

1=1 >=1 

As evident from our results in figure (5.4a) that reflection is of the order of 1% -4% and our 
results are in excellent agreement with those of Berenger in figure (5,4b) 

Here we would like to highlight the reasons for adapting to Perfectly Matched Layers. Zhang 
et al. [2] found that the Fourier transform of the time domain results is very sensitive to 
numerical errors, notably those resulting from imperfect treatment of the Absorbing boundary 
conditions used to truncate the numerical computation of an open structure. Thus, even 
though the time domain results may be reasonably accurate. The frequency -domain results 
obtained from their Fourier transform may not be acceptable as useful data. 

Even Zhang et al. [2] used a new super absorption technique for their research. This 
prompted us to look for other Radiation Boundary Conditions. We experimented with 
Perfectly Matched Layers and found its performance extremely superior compared to 
whatever RBCs discussed here. Details of PML and its Variant MIPML [32], using which we 
have developed simulator MIPML, have been taken up in next chapter. 
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Chapter 6 

Perfectly Matched Layers 


Introduction 

One of the basic considerations in any open boundary problem in FDTD is the reflection 
less truncation of the computing domain and also reducing size of the computing domain 
by placing it as close as possible to the scatterer under analysis. One of the earlier 
techniques [3], consisted of surrounding working domain with an absorbing material 
whose impedance was matched to the working domain. But this had the basic 
assumption that the waves are incident normally upon the interface of the working 
domain and absorbing material for reflection less transmission. So the absorbing layer 
had to be place quite far from the scatterer to make outgoing waves nearly normal 
incident. Another technique, which was discussed in previous chapter and has been very 
popular so far, is approximation of one way wave equation initially exhibited for 
acoustic wave by Enquist and Mazda [38]. Mur [39] and various researchers with higher 
order approximations applied this one way wave equation for creating useful ABCs. 
Though accuracy was reasonably good in time domain but for better accuracy in 
frequency domain, results in time domain needed to be further improved for some 
applications. Here, in most of the ABCs based upon approximation of one way 
equations, absorbing boundary needed to be placed away from the object for normal 
incident. Also second or higher order absorbing boundary conditions, which account for 
oblique incident, do not work for inhomogeneous structures because these have been 
defined for uniform space. 

Perfectly Matched Layers [13] do over come all above stated problems and accuracy is 
better than any other ABC by -lOOdB. It is Virtually reflection free for all angles of 
incidence, polarization and over ultra-wide frequency range [26]. 



6.1 Theory of Perfectly Matched Layers 

Maxwell Equations for TE to z case for a medium as shown in figure 6.1 having 
permeability as p and permitivity as 8 are given by 



Figureb. 1 E field and propagation vector kina medium 



Then impedance of above medium becomes equal to the impedance of the medium with 
same permittivity and permeability but with no loss terms i.e. p and a both are equal to 
zero. Impedance can be determined by first solving above equations to arrive at wave 
equation, which gives value of propagation constant and then putting value of it in above 
equations to find ratio of Ex/Hz or Ey/Hz, which gives the value of impedance. 
Hence 
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(6.3) 


(6.4) 


And solution of above equation is of the form given by equation (6.5) below, in a medium 
shown above in figure 6.1 

E = (.Eo COS0 u^, - Eo sin0 ( 6 . 5 ) 


After substituting above in equation (6.1) we get 


^ dy dx 


Using 

M = kcos9,k2 = IcsinO 
We get 
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(6.6 a) 


(6.6 b) 



Substitute value of k in above equation 

^ 

— = + j cos) (6.7) 

Hence 

Z ^ ^{jo)fx + p)l{j COS 4- O’ ] (6.8) 

If equation (6.2) is put in above Z reduces to 
Z = yIJiJs (6.9) 

Which is for the lossless medium. 

This proves that above medium (cT,p,s,|i) is matched to the irmer medium (0,0, E,p). 

Now lets first derive snell’s law and reflection coefficient for the interface of these two 
mediums, shown in figure 6.2 



Figure 6.2 Propagation of wave from working medium to absorbing 
medium 


Matching tangential components at interface x=0 we can write 

E,cose, = £, cos«, (6.10) 



Define 


(6.11a) 

(6.11b) 



El cose j cos0, 


Above is possible when 

/'clsinO,- =/:lsin0^ =^2sin0, 
This gives the snell's law 
sin0j- _ k2 
sin0^ k\ 


( 6 . 12 ) 


For a matched medium defined by equation (6.2), above law becomes 
sin0/ 


sin0, 

Hence 
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1 + - 


J(0^ 


0 ,- ^ 0 , 


(6.13) 


(6.14) 


For a matched medium incident angle is not equal to transmitted angle. 
Also from equation (6.1 lb) 

(l + r)cos0,- =rcos0^ (6.15) 

Since Hz is parallel to Y-Z plane. 

^ zi ^ zr ~ ^ zt 

zl z2 

z\ zl 

4(>-r')=4^ 

zl zl 


From equations (6.15) and (6.16) 


r = 


ZlcosO, -Zlcos0,- 
Z2cos0, +Zlcos0,- 


(6.17) 



For a matched medium though Z1 may be equal to Z2 but 0, ^ 0t hence reflection 
coefficient is not equal to zero. 

For this reason as mentioned above, in PML, splitting of the fields is done. By splitting the 
field components reflection coefficient can be made zero along with impedance matching. In 
equation (6.1) Hz varies with x and y both. So it is sp lifted in two components, one of these 
vary with x-direction only and the other one varies in y-direction only. X-travelling waves 
are absorbed by the conductivities px (magnetic) and ax (electric). Similarly y travelling 
waves are absorbed by conductivity subscripted y i.e. the py and ay. 

If both couples of conductivity are matched then medium is reflection less if it surrounds 
a medium with same permeability and permittivity but zero conductivity and resistivity. 
Referring to figure (6.1), we can write Maxwell equations and various fields components 
in a medium with anisotropic conductivity, as following. 


dE^ ^ , O’ K^zx+-^z>.) 

dt ](08 oy 

(6. 1 8a) 


(6. 1 8b) 

dH„ „ SE 

dt ox 

(6. 1 8c) 

dE 

^ d; 

(6.l8d) 

E, = -E, sine 

(6.18e) 

E^ =-E^cosB 

(6.l8f) 


(6.18g) 


(6.18h) 



Write equation (6.18 a)-(6.18d) in frequency domain and solve them for a and (3, after 
putting values of fields from equation (6.18 e)-(6.18h) in them. We get following 
expressions 
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(6.19b) 

(6.19c) 

(6.19d) 

(6.19e) 


So a general equation for any field component becomes 
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(6.20) 


Above is the solution for propagating waves in PML and it is evident that decay along 
the direction of propagation is due to artificial conductivity introduced. 

6.2 Propagation between Two Perfectly Matched Layers 


Doing the analysis similar to done above we can find the snell’s law and reflection 
factor for the wave propagation between two perfectly matched mediums and these are 


as following 
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(6.21a) 

(6.21b) 

(6.21c) 



Where G2 and Gl are defined through equation (6. 19c) 

If 

CTv , = < t ,,2 ( 6 . 22 ) 

And both the mediums follow perfectly matched condition of equation (6.2) then 
For all frequencies and angle of incidents 

6i = Qt (6.23a) 

r = 0 (6.23b) 

6.3 Practical Conclusions of Theory of PML 

As a summary of this PML we can say that propagation between two PMLs, whose 
conductivities follow matching condition of equation (6.2), will be reflection less at 

1 . The interface normal to X-axis if Oyi =ay 2 

2. The interface normal to Y-axis if axi = 0 x 2 

3. Above still holds good if equal conductivities are zero , which is the case with 
vacuum and a perfectly matched layer e.g. vacuum is (0, 0,0,0) and medium is 
(0,0,ay,py) at the interface n,6rmal to Y-axis or medium 1 is (crx,px,0,0) and medium2 is 
(ax,px,ay,py) at the interface normal to Y-axis. 

6.4 PMLforFDTD 

Though theoretically there should be no reflection if above conditions are satisfied but 
unfortunately there is numerical reflection if conductivity is changed abruptly. So 
smooth transition of conductivity is done over several cells thus minimizing numerical 
reflections and last PML is terminated Avith PEC though not compulsory. In equation 
(6.20) above last two exponential are responsible for amplitude decay and one of them is 
unity so equation (6.20) reduces to following at distance d 

^COS^ 

y/{d) = \i/{Q^ , ( 6 . 24 ) 


It reflects back from perfectly conducting layer and travels a distance d again hence 
reflection coefficient is 
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(6.25) 



Since so far we assumed a constant with respect to d which is not the case with 
numerical computation where a varies with distance x. So above becomes 

(6.26) 

Where d is the total thickness of the PML. 

With selection of following profile for variation of a 

(6.27a) 

(6.27b) 

Now it is quite straightforward to discretize equations (6.18a-6.18d) 

6.5 Short Comings of PML 

At this stage it is worth considering equation for evanescent waves. Consider following 
diagram showing propagation of a wave in medium 2 
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Figure 6.3 Propagation of evanescent wave in medium2 


Medium 2 



If E is an evanescent wave then angle 9 is imaginary. Here in direction x propagation 
constant is real and in y direction it is imaginary. Field components in the form of matnx 
can be expressed as following 

\Ex H 2 y ] = [e^ cosh 0 jE^ sinh 0 '\^ja>ii-ax-p y) 28a ) 


Putting above equations inequation (6.18a-6.18d) 
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Solving above we get 


-coaEg sinh0 


a = - 




f 1 


/ 


— II j + ~ ]sinh0 
cG jv (OS J 


V 


[cG, 


\ 


1-7- 


0)S 


cosh0 


Where 

G^ = JVy cosh^ e-W^ sinh^ 0 


W. = 


cr^ + jcos 


W 

V 

£ 


f , 


s IPx+J(^pJ 

CTy + jcos 

^Py^j(OP J 


(6.28b) 
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When values of a and p are put in equation (6.28a), above gives the field variation for all 
components. It can be easily seen that equation (6.28b) which express propagation along 
X direction, evanescent waves have natural attenuation which does not contain ox and it 
is not accelerated in PML. 

6.6 Numerical Experiments 

We have repeated the experiment done by us in chapter 5 with Mur s second order ^ 
approximation of one way wave equation as Absorbing Boundary Coiidition but this time 
with split formulation of PML proposed by Berenger [13]. This Experiment has been 
picked up from Berenger’ s original paper [13] and our results using program PMLTE , 
developed by us, are in excellent confirmation with Berenger. We have also repeated 
experiment done by Katz et al. [25] on the same software and results are identical. 
Reflection coefficients of the order of— lOOdB or high is quite normal which is — 60dB 



better than what Mur’s ABC could achieve at its best. If Thickness of PML is increased 
with optimum profile then any value of reflection coefficient is possible. In the program 
Developed by us we have used exponential differencing scheme as [13] and [25] both 
claimed it is a necessity because of sharp variations in the fields in PML regions. But it is 
not necessary to use exponential differencing as it produces same results as central 
differencing scheme does [15]. In our all experiments in subsequent chapters we have 
used central differencing and found it very effective and attractive because it is simple 
compared to exponential differencing. Our results and Those of Berenger and Katz have 
been enclosed for comparison. 

Also it is very evident that this material absorbing boundary condition is consisting of 
non-Maxwellian material and requires splitting of the fields. It leads to more computer 
memory than expected firom any Maxwellian formulation. Also as most of the time 
people work with Maxwellian formulations of the fields it is a bit difficult to understand 
this different formulation. For these reasons we have opted for unsplit formulation. 
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Figure 6.4 Berenger’s result for PML (4,L, 1) 
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Figure6.5 results produced by program ‘PMLTE’ for the same 
PML parameters as in Figure6.4 
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Chapter 7 

Unsplit Formulation of Perfectly Matched Layers 


Introduction 

In previous chapter, we introduced concept of Perfectly Matched Layers as Absorbing 
boundary condition. Field equations were non-Maxwellian in the perfectly matched 
layers and fields were to be splited up. From practical considerations of computing, it 
amounted to more memory requirements. Secondly splitting of fields introduced a sort of 
uneasiness because of habit of dealing with Maxwellian formulations. Following the 
work of Berenger [13], Chew and Weedon [18] introduced a complex stretched - 
coordinate formulation (cscf) of the PML, which was more compact. The principal 
advantage of cscf is the ease of mathematically manipulating the PML equations and 
providing a pathway to map PML in other coordinate systems such as Cylindrical and 
spherical coordinates. 

Following Chew and Weedon’ s work. Sacks et al. [19] introduced an anisotropic 
perfectly matched layer formulation in which medium conductivities were tensors and no 
transformation was required. This medium followed more familiar Maxwell equations for 
fields and thus it was termed as UPML (Unsplit PML) as no splitting of the fields was 
required. 

Here we will give brief introduction of cscf, then we will present UPML and various 
forms of it to take care of lossy and dispersive mediums. In this chapter, we will confine 
to UPML for isotropic medium in rectangular coordinates, as this is the basis for our 
Material Independent Perfectly Matched Layers (MIPML). However, it is quite 
straightforward to extend the same to other coordinate systems. S.D. Gedney has given 
details of UPML. [20,21] 



7.1 A Stretched Coordinate Formulation 

Here to present Berenger’s equation [13] in non-split form, a coordinate mapping into 
complex space is introduced which is given as 

|sy(y'W ,z-^ js^Cz'Vz' (7.1) 

000 


Where it is assumed that PML parameters, Sj =l+o'/jco8, are continuous function along the 
axial directions. The partial derivatives in stretched coordinates are 
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Thus the V operator in complex space is defined as; 
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Maxwell’s equation in the complex coordinate stretched space 
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Above can be easily derived from equations (6.18a-6.18d). 


Inserting E from equation (7.5c) into equation (7.5d) will give following wave equation 
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Setting determinant of above matrix to zero 
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Upon solving, further 



Now let us consider fields within working volume with total wave (incident and reflected 
both) as following 
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Now consider PML medium with following transmitted magnetic field. 
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Setting 


Ely =E2y 


From equation (7.8) and (7.9), at interface, existing at x=0 gives 
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Putting x=0 and canceling common terms on both the sides, we get 



Where 


r = (l + r) 

Above equation is true only if following conditions are met 
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Consider equation (7.7), which is given below again 

If following condition is satisfied 
s~^=Sy=s, (7.11c) 

Then above equation (7.7) reduces to 

Which can be written as following, using equation (7.10) and dispersion relation for 
isotropic media 


Which is same as equation (7.1 lb) 
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the cx)ndition for reflection less transmission for TMzcase is same as equation (7.11c) 
If a plane wave is incident on a half space with interface at x=constant plane and 
composed of a uniaxial medium with permittivity and permeability tensors 
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Then the plane wave is purely transmitted into the uniaxial medium irrespective of 
polarization, frequency and angle of incident. 

In UPML,Ex is discontinuous across the x=constant boundary ,while Dx=eSx‘^Ex is 
continuous.In berengers medium both Ex and Dx are continuous.These two mediums 
have different gauss’s laws. 

In comer regions following similar kind of analysis the general tensor is determined to be 


= e 


0 0 
0 £2 0 

0 0 . £2 


SyS, 


* 

0 


M = M 


//, 0 0 

0 0 
0 0 //j 


5A 

o' 

0- 




y 

0 


s. 


s s 

X r 


y 

0 


5;,^x 


(7.13a) 


(7.13b) 



In working volume Sj =1 where i=x,y,z and in regions outside the working volume 5, >1 
whefe i is the axis along which wave is propagating . In this region parameters for other 
perpendicular directions are 5 perpendicular to i=l -In corner region all 5i>l. 

For inhomogeneous mediums also it can be demonstrated that UPML does its job. 

7.3 Theoretical Performance of The PML 

If boundary is assumed to be PEC wall, power reflects back into the interior FDTD 
region.The reflection coefficient for isotropic working medium is computed using a 
simple transmission line analysis and has an amplitude of 

= (7.14) 


Where 0 is the angle of incidence ,d is the thickness of the PML slab ,'n is the 
characteristic impedance of the reference material ,and a is the conductivity of the PML . 
In this discussion, we have assumed propagating waves but in case of evanscent waves 
there is no acceleration to attenifation. So if 
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Kx>l will increase decay rate of the evanescent waves. 


7.4 The Discrete Space 

In discrete space spatial scaling of conductivity profile for magnetic and electric fields is 
done. Usually most successful are polynomial and geometric scaling .polynomial scaling 
is simply 
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If real part of s is greater than one then it can also be spatially scaled using polynomial 
scalling.To this end we let 
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For geometric scaling ,the PML conductivity profile is defined as 
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Reflection coefficient for geometric scaling can be evaluated by integration over PML 
as for polynomial scaling and it is 
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similarly if K>l,it can be geometrically scaled as 
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Here Oo is determined from R(0) which is predetermined along with g and d. 
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7.5 Descritization Error 

Design of effective UPML requires balancing of theoretical reflection error and 
Discretization error. If Omax is small, then maximum reflection is due to its PEC backing. 
To avoid it and some times with the intention of keeping reflection coefficient small, Omax 
is chosen to be large. But because of FDTD approximation, then Discretization error is 



too large. Total reflection error could be several order of magnitudes higher than 
predicted theoretically by Omax and R (0) relation. 

The largest Discretization error occurs at interface. Thus, it is desirable to keep step 
discontinuity at interface to be kept to a minimum. For the same reason, spatial scaling of 
conductivity is done but if further rise of conductivity after the interface is too rapid then 
Reflections from layers near interface are significant. 

With extensive experimentation, it has been found that optimum value of Omax is 
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This expression has proven quite robust for many applications. Above value may prove 
quite large when PML terminates elongated structures or sources with large duration. For 
the latter case, large reflection errors can occur in late time. It can be proven with 
following 

From equation (6. 1 la) we get 



■(7.20a) 


In discrete space magnetic tensor and electric tensors are staggered by one half of a cell. 
Therefore at interface we have Sx* =1. Thus above equation reduces to 




By putting value of Sx in above we get 
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If (0 »Ox(0)/So, then F - 0. If © « Ox(0)/Eo then r=l. Therefore all low frequency 
components will reflect back. A cutoff frequency is defined 


(7.20c) 

(7.20d) 

From above two equations it is clear that time exposure of PML can be increased by 
decreasing CTx (0), which mean this, will reduce R (0). So optimal value of omax will be 
smaller than given by equation (7. 19) 

To circumvent this complex frequency shifted tensor is used as a choice for s tensor. This 
will be explained in coming sections. 


7.6 Efficient Implementation of UPML in FDTD 

For derivation of the finite difference time domain expressions, starting with the UPML 
equations, Ampere’s law is- expressed as 
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If transformation is done from frequency domain to time domain convolution of electric 
field with tensor will result which is computationally quite expensive. To circumvent it, 
following transformation is done. 
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Now putting values of s in the equation we get 
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Now above can be converted to standard FDTD equations. Consider transformation 
equation (7.2 lb 1) after substituting expressions for Si. 
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Multiplying both the sides by jco and transforming into the time domain leads to 
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Similar equations can be written for other field components. 


7.7 The Complex Frequency Shifted Tensor 

In equation (7.20) we demonstrated that one of the limitations of using (7.21a) is that, 
signals of long duration suffer from large reflections off the PML interface. This effect is 
circumvented by using the CFS tensor. We let 
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To derive time dependent formulation for ampere’s law from equations (7.21a) and 
equations (7.22), consider for example a x-normal boundary (sz=Sy=l). Now introducing 
following variables: 
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Then equation (7.21a) can be written in time dependent form 
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From equation (7.22) and equation (7.23a.) expressions for Ex and Ey can be written as 
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Similar expression can be derived for Ez. 

Putting values of s; from equation (7.22) into equation (7.21a) we can drive electric field 
Ex, Ey and Ez as mentioned below in the corner regions. The constitutive relations are 
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We see that this method requires the introduction of three auxiliary expressions and 
variables compared to one required by the previous method. This results in an increase in 
computer time and rnemory. So this formulation should be reserved for special cases. 

7.8 PML Termination for Conductive Media 

Ampere’s law in the conductive media 




^ are relative dielectric constant and ^ = /c, + — 

' JCOSo 



Similarly by replacing x with y and z, we can write equation for y and 2 component 
respectively. 



Chapter 8 

Material Independent Perfectly Matched Layers 
Introduction 

PML technique originally developed by Berenger [13] has been extended for absorbing 
Evanescent waves and lossy materials [16], A review of PML in chapter 6 with numerical 
experimentation and a review of UPML covering absorption of evanescent waves and 
lossy materials in chapter 7 have been presented. Obviously, to further enhance the 
capability of existing PML technique, extending it to cover anisotropic materials, merits 
attention. It was proven [31] that under certain circumstances no PML condition can be 
derived even for the simple case when only the arbitrary anisotropic dielectric is 
considered. However it was realized by us and Zhao [32] that above difficulties [31] can 
be overcome by using flux D and B instead of E and H. Use of flux instead of field 
quantities makes it material independent, and extends PML to arbitrary anisotropic 
material. Though the proof published by Zhao et al. [22] was given for arbitrary 
dielectric medium 2-D, the proof developed independently by us is for a material, which 
is arbitrary anisotropic in permeability and permittivity both. Our proof is applicable to 
split formulations in 2-D like [22]. For two-dimensional case our proof clearly 
demonstrates that there is no need to use D and B as said in [24]. Use of E and H equally 
works. In 3-D, use of D and B instead of E and H makes PML equations simple. It also 
makes matching conditions independent of permittivity and permeability tensors, so they 
are simplified and make algorithm independent of medium. In fact we have demonstrated 
ven clearly that UPML equations are inherently Material Independent if material 
p ' erties are not included in Maxwell equations by properly transforming fields and 
dating them indirectly through the transforming equations. 



8.1 PML in Two Dimensional Anisotropic Space 

Expressing Maxwell equations as 
jcoD-Vy-H (8.1a) 

(8.1b) 

(8.1c) 


-jcoB = VxE 
D = £~^E 


— 1 

B = p H 
H. 




jiPx^ +Pyy+fiz ^ ) 
rr ~j^x>:+Pyy+fiz^) 


(8. Id) 


^ OX^ 


ov 


^ oz^ 


(8,le) 


Putting above in Maxwell equations will yield 
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By assuming permeability and permittivity tensors as following 
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Taking a two dimensional X-Y plane at z=0 where 
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In working medium 
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In perfectly matched medium. 
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In above equations value of propagation constant matrix is taken from equations (7.4) 
which are for Berenger’s split formulation. Matching tangential y components gives: 
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By putting phase matching condition in above 
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Where 
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Here we see that wave attenuates in x direction in PML. Though the splitting of fields is 
not evident firom derivation but it comes into picture with equation (8.13) and expression 
for propagation constant matrix in equation (8.9). 

8.1 PML in 3-D for Anisotropic Medium 

Here we will discuss some fundamental issues related to wave propagation in anisotropic 
media in 3-D space. In anisotropic medium Ex is dependent on Dx, Dy and Dz. It is the 
case with Ey and Ez. Similar statement can be made for magnetic field and flux 
components. This implies that in an anisotropic medium it is not possible to have TE or 
TM modes propagating. We have only Hybrid modes propagating. So any PML designed 
should be proved for Hybrid modes i.e. it should be able to absorb Hybrid modes without 
any reflection. All the derivations so far done have been for TE and TM modes in 
isotropic medium and curl equations in PML for isotropic medium are 

VxT ^-j(oix~sir (8.18a) 

Vxjr = -j(oe~sY (8.18b) 



Which strongly suggests that PML should be able to match any medium with the same 
■ matching conditions. But in anisotropic medium we do not have propagating TE and TM 
modes so it is not possible to prove UPML conditions separately for propagating TE and 
TM modes, and then assume that it will hold good for all propagating hybrid modes also. 
This has been the practice for isotropic medium but can not be used for anisotropic 
medium. We have tried to prove that UPML absorbs hybrid waves incident upon from 
anisotropic medium by numerical simulation. Prior to that we will demonstrate that 
UPML is essentially Material Independent and with the use of proper transformations it 
can in fact match any medium. Let us assume that Maxwell’s curl Equations, in PML for 
anisotropic medium, are 

VxE = -s— (8.19a) 

dt 

VxH = -s— (8.19b) 

dt 


Here 


^.x 


S = 


0 


0 


0 




0 


^2 


— — 
jO) 

y 

Sy =K+ — 
JO) 


s. 


<y, 

/jf -j ^ 


jo) 


( 8 . 20 ) 


ic>\ 



Kappa can be taken greater than one for absorbing evanescent waves. Also if medium is 
lossy, it can be taken greater than one to match it [21]. Its optimum value can be 
determined experimentally. 

From equation (8. 1 9b) and (8.20) 
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Before proceeding further we would like to bring out the special features of UPML 
equations. From equation (7.21c) after discretising we can write [20] 
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Referring to equation (7.17b) and taking 0=0 and Sr=l in this, we can write 
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If equation (8.23b) is put in expressions for A, B, C, D and E they become independent of 
So- So we can say A, B, C, D and E are independent of medium properties. If relative 
effective dielectric constant is not equal to one then it can be taken care by removing it 
and increasing value of m, reflection factor and d by the same factor. Which means that if 
any equation is of the same form that of equation (7.5) it can be matched by UPML. Here 
equations (8.19) are of the same form that of equation (7.5) with permitivity and 
permeability included in B and D. 

8.3 FDTD Equations for MIPML 

From equation (8.16d), amplitude of a wave can be written as below in a generalized 
form 

_ —COS e X 

\l/{x)=y/(0^ '' (8.24) 

Where interface of anisotropic medium and PML has been taken as reference. As is the 
usual practice either a PEC or PMC wall terminates a PML. Above mentioned wave 
travels from interface to PEC and reflects back and enters anisotropic media. So for a 
layer of thickness of d reflection factor can be defined as 



j^(^)^g(-2.Cr.c/cOS0/c) 


(8.25a) 


Agsin to minirnizs numerical reflections conductivity is gradually increased, with profile 
defined as below 







(8.25b) 


And Reflection factor is given as 


R{e)= 



(8.25c) 


Put equation (8.25b) for conductivity profile in equation (8.25c) and reflection factor 
becomes 
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(8.25d) 


Conductivity is given by 
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Here v has been replace by c. Though v is not equal to c but by appropriately choosing m 
and d the ratio between v and c is taken care of 
Putting in equation (8.21), following transformation 

SxD^b = s^Dx (8.27a) 

sPyb = SxDy (8.27b) 

s,D,b = SyD, (8.27c) 



We gel following equation 
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Discretising equation (8.27) and equation (8.28) we get 
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As discussed earlier we see that equations are identical to equations 8.21b and 8.21c. 
Also A, B, C, D and E are same as they were for UPML. 

Similar equations can be derived for Dy, Dz, Bx, By and Bz- Here first components, of B 
are calculated using the equations similar to equation (8.28) and equation (8.29). Then 
from B, Magnetic field H is calculated. From this value of H, flux density D is calculated 
which gives the E. This value of E is used to calculate B. Whole sequence is repeated for 
required number of time steps. 



8.4 Calculation of Field Components in MIPML 

In MIPML we calculate field components from flux densities. This is calculated as per 
following equations. 
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(8.30) 


8.5 Extension to Arbitrary Anisotropic Dielectric Cases 

The basic 3-D FDTD cell for MIPML is shown in Figure (8.1). It is inherent in FDTD 
that Dy and Dz do not exist at the same point where Dx exists. Therefore in equation 
(8.30) Dy and Dz are to be calculated at point (i+l/2,j, k). This is done by averaging Dy 
and Dz around the point (i+l/2J, k). Here averaging does not cause much of error as it is 
also second order accurate [20]. The averaging referred here is done in space on D and 
equation (8.30) can be written as following 
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By inspection equations can be written for remaining 5 field components. 
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Figure 8. 1 Basic MIPML CELL 

8.6 T reatment of Perfect Electric Conductors 

On the surface of PECs, the tangential electric field, E, must be zero, and the 
implementation of such a condition is straightforward in the FDTD updating formulation. 
However, from equation (8.21a) one can easily see that in the MIPML algorithm the PEC 
conditions Et=0 cannot be kept because E is updated from D. Therefore in order to make 
sure that the PEC condition is correctly imposed for Et, the values of the tangential 
electric displacement Dt on the surface of the PEC have to be calculated or justified. In 
following discussion. Different calculation procedures for Dt are adapted for the cases, 
first, when PEC is inside the mesh domain and, secondly, when it is at the boundary of 
the domain and exposed to incident waves coming directly on to it with no PML present. 
Both of these cases are present in microstrip problems. Here it is to be noted when waves 
come to PEC after passing through PML then D can be taken, as zero as it has been 



attenuated to a very small value and such an assumption does not introduce any 
significant error. 

8.6.1 If PEC is Located Inside The Mesh Domain 

If such a PEC is located inside the mesh domain, say at k = ko, then Dx and Dy can be 
calculated by imposing following conditions in equation (8.30) and similar equations. 
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Similar equation can be written for Dy. 

8.6.2 PEC at Domain Truncation 

At the domain truncation (k = 0), equation (8.33) can not be used as many of the points 
lie outside the domain. Using equations given below, we can first calculate Dx, Dy and Dz 
flux at k = 0 and put in equation (8.33). Though above equation requires only Dy and Dz, 
Dx is required in the corresponding equation for Dy. Equations for Dx, Dy and Dz fields on 
k=0 plane are 
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Put above equations in equation (8,33) to get following equation whose each point is 
known to us 

})’’ (/ + 0.5, ./,o) = - ^ * (2/);: (/, ./• - 0.5,1 ) - d; (/, j - 0.5,2) 

+ 2/j>;(/,./+o.5,i)-/2;;(/,j+o.5,2) 

+ 2/>;; (/• + 1, j - 0.5,1) - d ; q +\ j - 0.5,2) 

+ ID'l, {i + 1 , ,/■ + 0.5,2) - D; {i + 1, y + 0.5,2)) 

+ (/ + 1, ./,0. 5) - D" (/■ + l,y,l .5)) (8.35) 

Similar equation can be written for Dy. While writing code for any of these PECs in any 
plane it is always suggested to make space diagram to avoid any mistake. One sample 
diagram is given below. Similar diagrams can be constructed for other components in 
various PECs located in various domain-truncating planes. 



Figure8.2 components required for calculating Dx on PEC z 0 


Though condition of tangential electric components equal to zero at PEC is implicitly 
imposed through above formulation but it is always recommended that it should be 
applied on PEC wherever required. 


8.7 Numerical Experiment 

To validate numerically theory of UPML for arbitrary anisotropic medium we have done 
two experiments. 

For first experiment we will compare our data with Zhao[23]. For second experiment Mr. 
Zhao, Nokia Research center, Helsinki, Finland, provided us data from his paper “An 
Efficient FDl’D Algorithm for the analysis of Microstrip Patch antennas printed on a 
general anisotropic substrate”, which is written using split formulation. 


Both reference results have been obtained from split formulations. Our results have been 
obtained with unsplit formulation. In first experiment we have taken a domain uniformly 
filled with arbitrary anisotropic material with following specifications. 


Sr' No’" 

Parameter 

Value 

1. 

Number of cells (including PML)in 


1.1 

X-direction 

40 

1.2 

Y-direction 

40 

1.3 

Z-direction 

41 

2. 

Cell in PML in each direction 

8 

3. 

Number of time steps 

51 

4. 

Step size 


4.1 

X-direction 

1.3 CMS- 

4.2 

Y-direction 

1.3 CMS 

4.3 

Z-direction 

1.3 CMS 

5. 

Time step 

25 e-12 sec 

6. 

parallel Axis permittivity 

3 

7. 

Perpendicular axis permittivity 

2 

8. 

parallel Axis permeability 

3 

9. 

Perpendicular axis permeability 

2 

10. 

Optical axis angle with crystal X-axis 

30 

11. 

Reflection factor 

0.0001 

12. 

Conductivity profile index 

3 


— 

13. 

Number of cells in reference domain, in 


13.1 

1 X-direction 

110 

13.2 

! Y-direction 

110 

13.3 

! Z-direction 

111 

14. 

I vpc of excitation(Dz) 

1 

compact pulse 

15. 

1 Point of excitation in both of the domains 

i 

center 


In figure (8.3) reflection factor has been plotted. Figure (8.4) shows the results of [23]. 
Curve 2 is to be compared with our results. We find our results are excellent and 
comparable to those with split formulation. Here we would like to mention that due to 
limitations of computer memory our computing domain data is different from what has 
been taken in [23]. 

Reflection factor has been plotted on the interface (1 1:31,11,20) with MIPML, parallel to 
X-axis. Our MlPMl . performs better than -75 dB. 

In second experiment we work out characteristics of a Patch antenna printed on a 
anisotropic dielectric substrate with material parameter as given below in table. This 
problem has been taken for two reasons, first to demonstrate that our simulator is capable 
of taking up any planar geometry and secondly results for this geometry are available to 
us. 

Figure (8.5) shows the patch antenna geometry. Using the technique discussed above, this 
line-fed microstrip patch antenna printed on a general uniaxial anisotropic dielectric 
substrate is studied. The optical axis lies in the xz-plane. Following table gives 
description of the problem. 



r. No 

Parameter — 

Value 

i 

1 

Number ol cells (including PML)in 


i 

.1 ' 

X-direction 

58 

.2 

Y-dircction 

80 

.3 

/.-direction 

22 

: 

Number of cells in PML in 


:.l 

X-direction[ PML. left of the interface PML right of the interface] 

8 8 

:.2 

Y-direction[ PML left of the interface PML right of the interface] 

8 8 

:.3 

Z-direction{ PML left of the interface PML right of the interface] 

0 8 

i. 

Number of time steps 

8000 

1. 

Step size 


1.1 

X-direction 

3.891e-004 

1.2 

Y-direction 

0.4e-004 

1.3 

/-direction 

0.1985e-004 

). 

j 

Time step 

0.4e-12 

3 . 

parallel Axis permittivity 

2.31 

7. 

Perpendicular axis permittivity 

2.19 

?. 

parallel Axis permeability 

1 

?. 

Perpendicular axis permeability 

1 

10. 

Optical axis angle with crystal X-axis 

45“ 

11. 

Reflection factor 

0.001% 

12. 

Conductivity profile index 

4 

13. 

Excitation at two cells from interface with MIPML 

Gaussian pulse 

14. 

Halfwidth(T) 

15 picosecs. 

15. 

Time delay 

3.1 T 

16. 

Distance of reference port from patch edge 

10 cells 


rhe rectangular patch is 32Ax X 40Ay; the microstrip feeding port is 6 Ax X 27 Ay and 


jubstrate thickness is 4 Az. 


igure (8.6) shuws incident and reflected waveforms with respect to time. Incident 
,aveform is achieved by assuming no discontinuity i.e. no patch and reflected wave form 
1 achieved when patch antenna is present. Figure (8.7) shows the results supplied by Mr. 
hao. Figure (H 8) shows the Sn of the patch antenna and figure (8.9) shows the results 
f Mr. Zhao It seen that there is hardly any difference. This proves that simulator 
eveloped by us works correctly for arbitrary anisotropic media. 

1.8 C'onchisions 

code 'MlPMl .' has been developed for analyzing microwave circuits by extending 
nsplit formulatiims to arbitrary anisotropic materials. This simulator for designing and 
eveloping microwave circuits uses material absorbing boundary condition in which 
onductivities arc independent of material characteristics. We have clearly demonstrated 
hat all the paiamcters used in UPML formulation equations are indeed material 
adependent and UPML for isotropic materials can easily be used for all kind of 
aatcrials, Only parameters, Si, need to be defined properly. Further due to this special 
eature of MIPML, it can also be used to absorb waves propagating in material consisting 
(floss, dispersion, nonlinearly with slight modifications. We have also analyzed a planar 
liscontinuity of patch antenna type and found our results matching with those of other 
esearchers. This proves that software MIPML is correctly designed for analyzing any 
:ind of geometry on arbitrary anisotropic substrates and It can be used for the study of 
nicrostrip dispersion characteristics on anisotropic substrates with high accuracy. 
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Figure 8.3 Reflection factor by Simulator MIPML 
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Figure 8 b Results obtained by Simulator MIPML for propagation of the pulses 

e incident pulse is the reference pulse through the microstrip. Reflected pulse is the 
'eform of the reflected field from patch antenna. 
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Chapter 9 

Microstrip Transmission Lines on Anisotropic Substrates^ 
MIPML Analysis 


Introduction 


DciTiand lor high-speed pulse circuits has increased need for highly accurate full wave 
analysis of microstrip transmission lines. Till 1978 analysis were mostly TEM in nature 
for anisotropic substrates [34, 35]. Disadvantage of TEM analysis is that it does not take 
dispersion into account. These effects have to be included empirically. Later Full wave 
analysis method based upon waveguide model appeared [6]. Disadvantage of this method 
was that it did not take into account radiation effect and the surface wave generation. 
Also mode-matching step involved error because modes generated in model are not the 
same those generated in microstrip and microstrip discontinuities. For Anisotropic 
substrates, in 1985, approach based on TLM method appeared [36]. Also Koike [37] 
using Bergeron’s method presented his analysis. TLM is not as efficient as FDTD method 
is and with the invention of PML, FDTD results are far more accurate than those 
presented in [35] and [36]. Disadvantages of These methods have been discussed in 
Chapter one. 

Salient features of our EM simulator are that it is using MIPML, so capable of analyzing 
any material [chapter 8], Virtually reflection free domain truncation [chapter 8] and this 
reflection free boundary has certain advantages, which have been highlighted in next 
section and any kind of planar circuit can be analyzed without any modification to basic 
algorithm. For a new geometry, its dimensions are to be specified in subroutine 
DATAIN, implement boundary conditions for D in Subroutine CONDBND, implement 
boundary condition for E in Subroutine EFIELDS. The same routine also has to be 
modified for specifying dielectric constant in computing domain so fields may be 
calculated from flux. In subroutine EXCITATION we have to specify distribution of 
stimulus for the given geometry. We have already included a wide range of sources so 
any kind of stimulus just needs to be specified by entering choice. 
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For calculating dispersion characteristics of relative effective dielectric constant (REDC) 
of niicrostrip Zhang ct al. [42] have used multiple cell method (MCM). We also use the 
same and also proposed a new method of single cell (SCM), which is based upon 
transmission line analogy and contrary to the method proposed by liou [12], our method 
docs not require calculation ol'flux and charge. 

9.1 Methods of Calculating and Effective Dielectric 
Constant 


Koike [37] used following definition for characteristic impedance. 
Z„ = 2PJll (9.1) 


For evaluating REDC he used ratio of free space wavelength to guide wavelength. This 
was possible because his analysis was at single frequency. Though power definition can 
be used for signals containing wide frequency range. The only problem is that Ex, Ey, H x 
and I ly over a complete xy plane are to be stored and poynting vector is to be calculated 
by integrating sum of Ex*Hy and Ey*(-Hx) over xy plane and averaging over one period is 
to be done. 

Definition used by Zhang [2] for characteristic impedance is 

7 -EM 


(9.2) 


And for REDC, they took the Fourier transforms of Ex (t) at two different locations Qust 
underneath the center of strip), with a separation of L, along the propagation direction: 
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(9.3b) 


Taking ratio of (9.3a) and (9.3b), we get transfer function of this section of transmission 
line, which is 
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(9.4a) 


y{o))=a{(o)+ jf^{(o) ■ (9 4b) 

fi{(t)) = (9.4c) 

SrcJ! = (9-4d 

Results of both approaches in FDTD method are almost identical but second approach is 
simpler. It is because evaluating V and I is much simpler. This method is multiple z 
method as propagation of the wave is to be monitored at many locations along direction 
of propagation (z-direction). Advantage of this method is that results are accurate over a 
large frequency band, as it does not rely upon TEM approximation. Therefore, for 
frequencies where deviation from equivalent circuit model, from transmission line 
analogy, takes place results provided by the method used by Zhang et al. [42] are still 
acceptable. 

9.2 Shortcomings of Above Methods and Single Cell Method 

Biggest disadvantage in FDTD analysis is that the frequency dependence (dispersion 
relation) of the microwave parameters of the wave characteristic impedance and the 
relative effective dielectric constant (REDC) show oscillatory behavior. This behavior 
may be attributed to following reasons: 

(A) Imperfect ABCs that generate artificial reflection at computational boimdary 

(B) The excitation source can contain mixed modes of propagation and evanescent waves 

(C) The continuous radiation loss spectrum that was taken into account in FDTD, but 
failed to be implemented properly in the multiple-z method for phase velocity 
determination 

(D) Possible numerical errors due to the numerical approximation of the Fourier 
integration. 

The first reason is supposed to be the biggest contributor to the oscillatory behavior. By 
using MIPML formulation we have reduced it. Our results for multiple z method show no 
oscillatory behavior. A slight droop at very low frequency is observed in REDC. We used 
a modified single cell model to correct this small droop to get perfect charactenstics. The 
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other biggest nuisance in multiple z method is due to the multiple of an integer with In in 
phase ditterence, when trying to determine a unique propagation constant. Liou et al. [12] 
have tried to overcome above problem by introducing a new method. 

I he method used by Liou [12] requires charge Q and magnetic flux for evaluation of L 
and C' per unit length of transmission line. He used following fonnulation 
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We use different method of evaluating L and C. we have straightaway discretised 
transmission line equations. 
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Once L and C are known we use equations (9.5d) and (9.5e) for REDO and Zo. We have 
observed distinctive improvement in REDO results. Zo results are same as those worked 
out with V/I method. It shows accuracy of this method . Here caution is to be exercised 
that for frequencies where equivalent circuit approach holds, results obtained by a single 
cell method arc accurate but for higher frequencies, where assumptions of equivalent 
circuit approach deviate multiple z method gives accurate results. For this reason we have 

preffered to give results by both method. 


9.3 Results 


Wc hti'vc enclosed results troni Rcinmut K. Hoffnian for Crystal iriain axis perpendicular 
to the substrate surface. These results have been calculated by the methods of lines with 
equivalent transformation in figure (9.2). In the same figure we have compared our 
results. We have enclosed our results for various w/h ratios. It is observed that results for 
characteristic impedance arc in excellent agreement with published results. There is slight 
discrepancy in RTDC but this is owed to the fact that our results are for full wave 
analysis and those given by I loffman are calculated with TEM approximation. To rule 
out possibility of spurious modes affecting our results, which have been calculated by 
sampling fields at approximately 1 .5 A., we compare them with results calculated at 4 k. 


Sr. No. 

Method 

Zo 

REDC 

1. 

I loffman’s results 

46.5000 

7.3750 

2. 

Results sampled at 1 .5 A 



2.1 

j 

MCM method 

43.0943 

7.7240 

2.2 

SCM method 

43.4312 

7.6734 

3. 

Results sampled at 4 A 



3.1 

MCM method 

43.1011 ! 

7.5487 

3.2 

SCM method 

43.1028 

7.5451 


In above table we see that results at 4 A, by MCM and SCM are exactly same. The results 
are not very different from the results at 1.5 A. So results at 1.5 A are , well within 
reasonable accuracy. This reduction of computing domain enables MIPML to be run on a 
Personal computer and still get accurate results. We have also worked out effect of angle 
of optical axis with respect to X-axis in figure (9.3). Results of Koike [37] have also been 
enclosed figure (9.4). Our results show the same pattern of variation. Details of structure 
are shown in figure (9.1). SIMULATOR MIPML produces results, which are most 
accurate of all. This is due to most accurate boundary conditions and because of the use 
of single cell method. These reasons have already been mentioned above in details. We 
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have enclosed results at 1 .5 X though the results at 4 X or farther would have been more 
accurate. I he reason is that tor ver}’ odd w/h ratios the size of computational domain 
would have become so large that it would have been impossible for us to calculate results 
on the compute I available to us. lo maintain uniformity in the results for all w/h ratios 
we have calculated them at 1.5 X. I*.vcn there also our results are very accurate. For any 
specific applications results can be calculated with larger domain. 

It is demonstrated in the results for w/h ratio of 0.8 (figure 9.6) that nuisance is faced in 
MC'M because ol' multiplication of an integer with 2?! in the phase difference when trying 
to determine a unique phase difference. So selection of sampling points should be 
carefully done. In SCM, We demonstrate in the results mentioned above that there is no 
such problem. This gives capability to calculate results upto a larger bandwidth with the 
same cell size and time step. This is quite a significant improvement. The same is again 
demonstrated in figure (9.7) that if sampling points are not properly selected in MCM 
than results could be disastrous. But SCM does not pose this problem. From figures 
(9.1 1 ) and (9.12) capability of SCM to reduce oscillatory behavior is very clear. These 
properties of SCM make it quite attractive. 

9.4 Conclusions 

Though we have tried to compute results for a wide range of w/h ratio which has never 
been attempted earlier. The reason is for very odd ratio requirement of a large computing 
domain makes it very difficult to do it with a reasonable amount of computer memory. 
Also increased computing time compounds the problem. These difficulties come into 
picture because of adequate sampling of microstrip geometry dimensions make it 
necessary to keep size of space sample very small. Also sampling needs to be done atleast 
one wavelength away for spurious modes to die down. These are generated because it is 
not always possible to give excitation in pure dominant mode and because of 
inhomogeneous medium, evanescent waves are generated. If sampling points are kept at 
one wavelength away with very small space step due to geometrical constraints, then we 
have a large computing domain. Also thickness of substrate or width of strip should be a 
fraction of wavelength and this requirement also leads to above mentioned problem of 
taking very small steps which may not be required for the given frequency range and this 
means a burden on computing resources. In fact it required extensive experimentation for 
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a very long period of time h\ us to determine optimum sizes of computing domain for 
various w/h ratios for the gi\ en computer memory. 



Figtme 9. 1 Gross section of rnioiostrip line along vtithexientatian of 
oyital axes vwtii respect to XZaxes. 


>1N 


Hoffman s results for dielectric constant 
Hoffman's results for 
Our results for using SCM 
Our results for Z^ with MCM 

Our rsults with MCM for REDC 
Our results for REDC with SCM 

Figure 9.2 


0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 3.6 4.0 4.4 4.8 


w/h ratio 









0 30 60 90 

0 (dejree) 



111 



























w/h=2.0 theta=0° 


40 


10.0 


38 


36 


° 34 


N 


32 


30 


28 


0 


9.0 


8.0 


5 10 15 

Frequency in GHz 


20 



Figure 9.9 


116 


REDC 




































Conclusions and Further Developments 


In this thesis work we have proved that any field quantities, if they can be expressed in 
Maxwellian form, can be matched using unsplit formulation of PML. This is possible due to 
material independent nature of the PML equations. We have also demonstrated that in 2-D space 
any incident plane wave irrespective of its polarization, frequency and angle of incidence can be 
matched, at the interface of PML and Anisotropic medium without any reflection. Our results 
clearly demonstrated that because of virtually transparent boundary conditions, there is no 
oscillatory behavior in frequency dependent characteristics of transmission line and results are 
highly accurate. Since by using MIPML formulation it is possible to truncate computing domain 
(compare the size of domain used by Zhang et al. [2] and by us for similar w/h ratio), by using 
reasonable sizes of computing domains, very small w/h or very large w/h ratios can be analyzed. 
For example we have analyzed w/h ratios from 0.2 to 5. We could have gone for even smaller or 
bigger ratios then this if computer system with slightly larger memory were available. 

After formulating MIPML for anisotropic media and preparing simulator for it, we tested it 
thoroughly for various planar circuits and for a uniform 3-D anisotropic space. We could achieve 
very good absorption with worst case figure well above -75 dB in case of uniform arbitrary 
anisotropic media. 

We also simulated propagation of a Gaussian pulse and calculated characteristics of a patch 
antenna. Our results were in excellent agreement with Zhao who calculated them using split 
formulation. 

After making sure that simulator MIPML properly and gives accurate characteristics we used it 
for our ultimate goal of finding dispersion characteristics of Microstrip on an anisotropic 
substrate. We have chosen it for several reasons. One of the main reasons was that most of the 
researchers [34, 35] so far used TEM mode approximation as the mode of propagation while 
calculating characteristic impedance. They did not give dispersion characteristics in their results. 
As frequency increases mode of propagation does not remain TEM. Mariki [36] calculated 
dispersion of characteristic impedance for w/h ratios of 1. 5,3,3. 5 and 5 but as far as relative 
effective dielectric constant was concerned he gave characteristics only for single ratio of one. 
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We have worked out ratios as wide as 0.2 to 5 and for all these ratios, we have given dispersion 
characteristics of relative dielectric constant and characteristics impedance. This data is very 
useful for designing microstrip circuits on anisotropic substrate. 

We have also worked out effect of angle of crystal optical axis with X-axis. Thus effect of angle 
has also been taken into account. This shows that our simulator is capable of working and 
providing accurate results for microwave CAD in non-homogeneous, arbitrary anisotropic 
media. 

We have also used a new method, which is very simple and provides highly accurate data, for 
dispersion characteristics by reducing oscillatory behavior. Here we use data over a single cell. 
So nature of errors in time domain which are already very little because of excellent boundary 
conditions is identical and magnitude is almost same as they are over a cell only. Since we are 
Subtracting voltages and currents, these errors also tend to get cancelled. This accounts for non- 
oscillatory behavior. Our method is simpler than liou et al. [12] as it does not require calculation 
of charge and flux. Also conceptually it is very simple, as it is just descretization of equivalent 
circuit of transmission line analogy. However, pinch of caution is that in frequency range, as 
long as above-mentioned model is valid, results of this method are valid. For further frequency 
range above 25 GHz, multiple cell method [2] is to be used. 

Further Developments and Suggestions 

In this simulator we have kept provision of taking care of evanescent waves and lossy media. To 
do so Kappa factor in DAT AIN subroutine has to be kept more than one. Though we did not do 
any experimentation with this because of shortage of time but it can further reduce size of 
computing domain, as we need not to provide space for attenuation of evanescent waves- Some 
very exciting experimentation will be done in near fiiture for calculating S parameters of basic 
microwave discontinuities on arbitrary anisotropic substrate with very small computing domain. 
As we have already mentioned that this simulator can analyze any planar geometry without any 
modification. 

We also plan to analyze scatterers, which are not planar. For that slight modifications in the 
subroutines responsible for implementing boundary conditions for D and E are to be done- There 
is no need for doing any other modification in the core algorithm. With this addition it will be 
possible, to analyze DR circuits. Ferrite based circulator etc. 
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most exciting area could be the analysis of multiconductor microstrip and stripline circuits, 
hot topic in the area of computer hardware. Our simulator can do that without any 
edification. Only geometry needs to be specified. 

I research the trend is to build upon what others have done. Though our simulator has been 
sted in Cartesian Coordinates but as we stated earlier that shortage of computer memory, 
)WSoever big computer may be, always comes into picture. In our case we faced it for very 
nail or very large w/h ratios. It will become a very potent simulator if in future, Subcellular, 
ubgridding techniques [4,5] and conformal coordinate [6] techniques in form of new 
ibroutines are included with this simulator. This will facilitate doing simulation for thin gaps, 
lints and wires. Also cylindrical DRs and cylindrical ferrite geometry based circulators can be 
esigned. This simulator, named MIPML will become a complete package for any kind of 
eometry 

astly we would like to stress that this simulator does simulation for isotropic media also. This 
mulator can also do any problem of any planar circuit in any linear media including anisotropic 
ledium. 

ince this simulator has been developed using Material Independent properties of PML we have 
amed it MIPML Simulator. 
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